{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the basic libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython import display\n",
    "import os\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the environment for running the experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paths for loading/storing results.\n",
    "name = 'results/Gaussian_MINEE' # filename\n",
    "chkpt_name = name+'.pt'              # checkpoint\n",
    "fig_name = name+'.pdf'               # output figure\n",
    "\n",
    "# use GPU if available\n",
    "if torch.cuda.is_available(): \n",
    "    torch.set_default_tensor_type(torch.cuda.FloatTensor)\n",
    "else:\n",
    "    torch.set_default_tensor_type(torch.FloatTensor)\n",
    "\n",
    "# initialize random seed\n",
    "np.random.seed(0)\n",
    "torch.manual_seed(0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data using the gaussian model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from data.gaussian import Gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_size = 1000   # sample size\n",
    "rho = 0.9           # model parameter\n",
    "\n",
    "rep = 1             # number of repeated runs\n",
    "d = 6               # number of dimensions for X (and Y)\n",
    "\n",
    "X = np.zeros((rep,sample_size,d))\n",
    "Y = np.zeros((rep,sample_size,d))\n",
    "g = Gaussian(sample_size=sample_size,rho=rho)\n",
    "for i in range(rep):\n",
    "    for j in range(d):\n",
    "        data = g.data\n",
    "        X[i,:,j] = data[:,0]\n",
    "        Y[i,:,j] = data[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A plot of the first dimension of $Y$ against that of $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2df7SdZXXnvztcco2BVG9hGlHujY2U1sm6S20GZY1VplJDCtWW1o6ty9hx1Szq2KldmjglKCqhC8gsll214zIzUon1FzPE6iLQoONQxmmwJi6MQXCGWG5ACKDXNBAzide7549z9r37PPd5f53znvO+57zfz1pZ3HPP+77P8wbY+3n23s93i6qCEEJI81hW9QQIIYRUAx0AIYQ0FDoAQghpKHQAhBDSUOgACCGkodABEEJIQ6EDIANHRO4RkT8c0Fh/JCJPisizIvKzOa5/REQuHcTc6oCIfFJEtlc9D1INdACkL7QN6cm24X1SRP5aRM4q+Iw1IqIiMtblHM4EcDOA16vqWar6w26ek/J8FZGXlPlMQgYJHQDpJ7+hqmcBeAWAfwXgmgGP/3MAngPggQGPS8hQQAdA+o6qfh/AXQDWhd+JyDIRuUZEZkTkKRHZJSI/0/763vY/j7V3EhdH7h8XkY+IyOPtPx9p/+4XAHzX3f/V2NxE5K3tsX8oItuC7y4SkX0ickxEnhCRj4rI8vZ3Nrdvtef2b0Xk+SJyh4g8LSI/av/8oqS/FxF5n4h8X0SeEZHvisjrssZtf68i8k4R+b/te68TkbXte46LyG1unpeIyGMicrWI/KC9M3tLypyuEJH722P/g4hMZ82XDDGqyj/8U/ofAI8AuLT98/lorcKva3++B8Aftn9+O4CHAfw8gLMA7AbwqfZ3awAogLGUcT4M4D4A/wLAuQD+wY2Tej+AlwJ4FsBrAIyjFS6ac/P+ZQCvAjDWftaDAN7t7lcAL3GffxbAbwN4LoCzAfw3AH+bMPaFAB4FcJ6b69oC434JwCoA/xLAKQD/o/13+DMAvgPgbe1rL2m/083td3wtgBMALmx//0kA29s/vwLAUwBeCeAMAG9r/3scT5sv/wzvH+4ASD/5WxE5BuBrAP4ewJ9HrnkLgJtV9Xuq+iyAPwPw5gJx/7cA+LCqPqWqTwP4EIC35rz3dwDcoar3quopAO8HMG9fquoBVb1PVedU9REAH0fLgEZR1R+q6u2q+mNVfQbA9SnX/xQtw/pSETlTVR9R1cMFxr1RVY+r6gMADgG4u/13+M9o7bZeHlz/flU9pap/D2APgN+NzOkdAD6uql9X1Z+q6q1oOZdXpc2XDC90AKSf/KaqPk9Vp1T1nap6MnLNeQBm3OcZtFa+P5dzjNj95xW491H7oKonACwkikXkF9phnKMichwtB3ZO0sNE5Lki8vF2SOk4WiGs54nIGeG1qvowgHcD+CCAp0TkcyJyXoFxn3Q/n4x89gn3H7XfzUj6O5oC8J52+OdY23mfj9aqP3G+ZHihAyBV8zhahseYRCtk8SRaoY5u7n8859hPoGXgALQMOFphHONjAB4CcIGqrgJwNQBJed570AqVvLJ9/Wvs0bGLVfUzqvrq9vwVwI1djpvF80Vkpfuc9Hf0KIDr207b/jxXVT+bMV8ypNABkKr5LIA/FZEXt8tE/xzA51V1DsDTaIVkfj7j/mtE5FwROQfABwD8Tc6x/zuAK0Tk1e2k6YfR+f/E2QCOA3hWRH4RwB8F9z8ZzO1stFbfx0RkAsC1SQOLyIUi8qsiMg7g/7Xv+2nOcbvhQyKyXER+BcAVaOUnQv4LgKtE5JXSYqWIXC4iZ2fMlwwpdACkam4B8Cm0wiX/hJZx+WMAUNUfoxVH/9/tkMSrIvdvB7AfwEEA3wbwzfbvMmnHz/89gM+gtRv4EYDH3CXvBfD7AJ5Byzh+PnjEBwHc2p7b7wL4CIAVAH6AVmL671KGHwdwQ/vao2glsa/OOW5RjqL1bo8D+DSAq1T1ofAiVd2PVh7go+3rHwbwBznmS4YUUWVDGEJGFRG5BMDfqGpiOSppLtwBEEJIQ6EDIISQhsIQECGENBTuAAghpKF0pbJYFeecc46uWbOm6mkQQshQceDAgR+o6rnh74fKAaxZswb79++vehqEEDJUiMhM7PcMARFCSEOhAyCEkIZCB0AIIQ2FDoAQQhoKHQAhhDQUOgBCCGkodACEEDIAtuzahy279lU9jQ7oAAghpKEM1UEwQggZNmzVf3BmtuPzjk0XVzYno7IdgIg8R0T+UUS+JSIPiMiHqpoLIYQ0kSp3AKcA/KqqPisiZwL4mojcpar3VTgnQggpFVvp12nlb1TmALSlQ/1s++OZ7T/UpiaE1J46GvNuqDQHICJnADgA4CUA/kpVvx65ZjOAzQAwOTk52AkSQkhJ1NFZ1KIhjIg8D8AXAPyxqh5Kum79+vVKNVBCSFWECd3pqQkA9TTuHhE5oKrrw9/XogxUVY8BuAfAZRVPhRBCGkNlISARORfAT1T1mIisAHApgBurmg8hhGRR54RuN1SZA3gBgFvbeYBlAG5T1TsqnA8hhDSKKquADgJ4eVXjE0JItwz7yt+oRQ6AEEKGiTrq+nQDHQAhhDQUagERQkhO6qzr0w3cARBCGkeREM6ohHticAdACCE5YRkoIYQMKUVCOKMW7olBB0AIIQUZFSdQCy2gvFALiBDSK1t27cPho8exdvWqXIZ8FFb+tdYCIoQQEqefSWiGgAghlTHI1XUY07ffZY09zCv/LOgACCGkhgwiCU0HQAgZON0YN39NN8Zw1Eo4y4AOgBBCasggHBYdACFk4BQxbuFu4cqb9uLEqbnc9yeNnUUTdgp0AISQxjIMRr6fc6MDIIRURh7jFtstdGu4u9lx9MNJ1MXx0AEQQhpHE2Qe8sCTwISQkSc0+CvHW2tfyyVMT00A6HQA/Vz52zxi4/aDpJPA3AEQQhrH2tWrOj43beVvcAdACGkM4aq+qtDPoMflDoAQ0heG2Yg2deVv0AEQQkqjbsnUcD51mVdd5kEHQAhJpGjDFJNZHsSc/MEwANi9dUNfxx1FKnMAInI+gF0AVgOYB7BTVf+iqvkQQnrj8NHjOHFqDgdnZivfCbDMMx9V7gDmALxHVb8pImcDOCAiX1bV71Q4J0II8hlQn0i1lb+XWu4XNq6t/HuRhYiR9ZxRciaVOQBVfQLAE+2fnxGRBwG8EAAdACFDhnXXqotxrLrKZ1ioRRmoiKwBcC+Adap6PPhuM4DNADA5OfnLMzMzA58fIU2lqAEt2m6xV7LGKzL/rENaVR3iKoPaloGKyFkAbgfw7tD4A4Cq7gSwE2idAxjw9AghOYjp9NSBYTDOVVLpDkBEzgRwB4C9qnpz1vU8CEZIPYnlAXpdIYeOJCbTUHSsoo1nun1G3ajdDkBEBMAnADyYx/gTQvpLN+Gb0BAfPrpkE59rnDRjT/pHlSGgfw3grQC+LSL3t393tareWeGcCCEpZK1+165e1XMO4PDR4x1NX2JjJzmOrHnnKQttUpP4KquAvgZAqhqfENLiypv24uTpOcy3o8EHZ2Zx5U17cxnxNK3+EFv5m2E/ODOLjdv3YMXysSXGflnJliHPzqSJVJ4EJoTUn25X0N3Gy80ZmWxzljPKmhfVP+PQARDSUMxY2up75fgYTp6ew7rJicIGsmhnr0NHZrFi+diCfIMd6jLCHUFeQimKWB+Ak6fnsGXXPjoB0AEQQnJQ9GBVXr2emBMC8lf1xE4jp91jOYoyGMZqoBA6AEIaSiip0E8xtZOnO+Ua5rVl9GPG2Fbwh48ez7VS97kFy18ArfcJHQSA2ugV1QE6AELIAlkNU/LWxifp9Vjox/DhGr+CL3qYbFA6RMBoCc3RARAywuSRc/ZiamXLOdvKP8SHeMIV/MbtewC0dgl5Vure2Zw8PRcVh+t3r99hhQ6AEAIAS4ywL8+MlXn6n5NWwyuWj3U4Aavu8Y4mDAPNB+IEg+gxUIRREpqjAyBkBKmDnLONaU7ADHuYiLXPK8eTzwPkyQP4+a8cH0t0HMNssMtmWdUTIIQssmXXvkrE1MxQWhXOvLZ+t3J8DNNTEx3O4uDM7JLQzPTUxMJ1YVjnxKk5rFg+hmWChWvM6Nvz7LoQcxp1EpgzwtDSMMIdACEjSJ4wRVLC1zh89HhiDD/Er+jDE7/GiuVjHUY/vG+ZLA3/AIvVQn5+Se8X7mZY7ZMOdwCE1ABb+Xujlbbq7cdOYcemi7F76wasHB/DyvExrF29quNQmK14p6cmFnYK9t3a1as6wi3+Z3vW7q0bloRk1q5ehWXSuuauay7H9NTEQtjHy0GYE6CkQ7lwB0DICJO3mbu/1lb9eZqupz0rdjgrtnL3p4Dt88nTcwvlon4nEdb6p7WDtGu48k+GDoCQGpC3siTJ4ALIdRI2DTOYsTBMOL4ZYtsJJJE2HwvrxMpQLS+wcryVO1ixfCx3kjpW3grQEcSgAyCkYcROx4YrcU+ecwIx6Ya0cI3fIaQRykuHuxEa996gAyCkRhSVXwbQkXBNknJOk2sOV+JAckK225O6aViOIO2wluU8wlV9jFGq0+83TAITUnOKJnyTNHY8XhcnvN5CLlba6ZO+Nh9LyFrS2nPlTXsXGrqcODW38Dl8J18mak1gijqWoo1hSCfcARAyBITCaOGqNkycWpgmpsrpdXmMrNO6IfZ7X86Z97Ru2IAmRvh+RVb1XPnnp9Km8EVhU3jSJMJ6+lioxNi4fc8Sg+qbqfgTsqH+f6zaxsYywlW+xfyB5ORzUvVQWMFjOw77nCYFnSZD4efkoSOoYVN4Qkg2fmXuD0SFidy01TTQaVRDAwx0SjCnNWv35JVXDp/RTbOXPHr/fj4AMiuUCB0AIbUkybDnDbMUMX5pz4xVDIWJ4NipXmBppU4S6yY7V+6xZHDW3Dw2H5aBZkMHQMiQsEyWGkcfIvKVO7ZzMEllHx4J4/dG2tmDWFI5a5cQzm96amJB9M3PJewd4DE5ijwN65N6EJBk6AAIqRnhQSYz7PNabDXrdw+9SCrHdH0sDOX7++Yp0fSOyZ5tz4r18o0Z8Sx9IjaAzw8dACEl0GuYoWgzljQDCXQmVcPksRnutKRynjEAdISFbEfhu375ew8dmY2GtEKBuNjfhU9YZ512ZmlofiqtAhKRWwBcAeApVV2XdT2rgEhdKdMBhMlMC52Eydkw/GPG0eZh5ZY+xh47gRt7rj0nZrRj0gwW1vFGOtw5+BCVXb9764YluQTvmOy7MBTk/07CKiCu+JdS1yqgTwL4KIBdFc+DkK7otT9saMjtOV4JM0Ys/BI7ALZi+ViHIY1JNafp68Sqi+a1Zej9WKFgm30XOq/YqjzmfMK/v6S8BemNSh2Aqt4rImuqnAMhVROLac9rcjes0OmE5ZyhtIM//JVEeDjr4MwslsnSGn1/PqAXaebYO/mWlFklpezxWw5V7wAyEZHNADYDwOTkZMWzIU0lyciknVDNY5h8SMbCJ2nGOik+DuQvEY2R1o0rVOFMWsnbIbOV42NR6Wi7N4myW1KSbGrvAFR1J4CdQCsHUPF0CCmNcCVvK34z6El17L59o4V+rAeu7SbC5u1Z8gu20rdx/YnhcKUfKwtdJovzPXl6ruOwWuwdQvJKPcR+z5V/99TeARBSJXlj/LFGJEUPIlmyNhZaOXRktiOUE2ugDnT20LWdwpZd+zK1d6zz15U37cUyQUcSNzyNnLVKDw92Fa1wSoKhnvKhAyCkItL62hrmDGx1bvhqobWrV3UIuwGdydLDR493hJZiUs+Hjswu0ROy69PCUlaeGVYaGXlj+7F7Sf+p1AGIyGcBXALgHBF5DMC1qvqJKudEiCctNJHUnCTUounFqFm1TSjg5sfxTiLNyBtZukFAp1JoUmze7xT8PIxYxVE3O4FeK61IMlVXAf1eleMTMkiSVvoWetm4fQ/WTU501MPPa2eCNhbK8bX/G67bAwAd1Tx58PX54Wo/3F34e+z5edpDkvrBf2OE5CAtxj+INoVhpU9SC8e0lb9/lpeB9k4llGoA0kNAnljOI9wFdJMHYIev/kEHQEifiTVlARYNpo+9Hz56HMtkMSnrT8L6Z1hi14y4ddTK0zzdV+v4HUd4+Mw+7966obDAms1n7epVS/r6kvpAB0BIm7wrTFvpJzU8ybo/S8xsXhdF0sxw+lg7gAWp6FhVTt4DWjaO/+ydQJh4TsPu847LN5wvQ5eHzqN82BOYkD6zY9PF2LHpYqwcH8PK8TGsm5zAusnFfrt3XXP5wrW7t27A9NTEQqgkXDWvHB/LlIk4eXpu4RRvFqEzioWPtuzah91bNyzMzXoE21xsx7Jieesswsbte3BwZnbBOZmjpAGvH9wBkMbTbZVJ0mnX2POT9H7mtVMHP9baMXaI6+TpOUxPTURDPXadOQD7p63w/XOS8gW+2iiGF4tL+x2pN3QApFF0k0hMq94pIzGZFRKKVf6E4ZsYdo93COsmJxbuiyV8jROn5hacRax+3z6H1Uahg7FSUa7+6wlDQKTxWIjGwhv2uQzSNPVtRb5i+dhClU8s7h6erDVOnm5JNux9/2IIyUIysfCPOQ2bizWZScJ/x1X9aMIdAGkE3YR5sqp3vDRD2mneNGJ1/iGHjsxG6/NNljk2lq/u8Vh+IO1wV8wphLLSvqlMklIoV/71hw6AjDxJPW1D+mGs8vap9W0fQ8xJJBnnsJIHSD4A5u+3OH+enIC9Q3iqNyt8ReoNHQBpBD5BaavZJNJi/uHnpO/C9oY2h7Dbl5G0YvckOQegM9Gbh1if3zRCQ2/vmaYBROoPHQAZWWJaNEDL+NohpW6SwVnEdhu+I5Yd9gpj7HlO8SaR9z7fzjHJCcQSueHz0yQgylL/JP2HSWDSSExkbcuufR3tEq05uVW+eF17f7AJaJVs2u9C52B18sDiYS27pkzDaPH8LB0e393LnNHurRs6avmNok7IDqoBcfVP7gbqC3cAZCRJUuYEkg1wUp7Aeu/6kkg7ietDI/46oGVUY89M64GbpPmftDuYVywZM2y8bljzGBvH/h7CE8B5iJ1BKEv9kwwOOgDSKEJ5ZSO2OvchG298veHzipihkbbVNrB09xBiOw+Pn2NeQTYv8AZ06vh4Tp6e63AcoZS0vVtRfHWRVQWxEqi+0AGQkSSrHDOmn2NhC1s524q+iCaO4Ru6J82h16SpORxvvC22752P4U/3xiSlfbLaOz9v1KenJpa8V1KDdxr/+kMHQEaCvKdzw7JMr1RpRs4bzlBd04dVwlBNbJWepyNWbIUektb+0RNKToSlm3aNx4eqyjLaNP7DAR0AaQThAbAYtroN7zHjuEw6Nf8Pzsx2dMXyjsIkF2I5AEs+l2Ug00JDSXX6PlxlO5zQaGc1dWdbx+GHDoAMNbHTuv50bnhYySQNisS3fdhk4/Y9WLF8DLu3bkhczdv1YbglbOpimFOxbl4x8paI5lH3NPx8fOUTQCPeFOgASGWUIcyWha3OrflJeBDLh2fC+LaFdsKTsmZkY3OwBu1h4tjG9LmH8JBYFnmcQF7HZge48vYOoEMYTegAyFCS53SuGf6i/XH9vTFMwtknS8PvfVjGG27/syWYrTKpbPa+//IlpalAPD6f17FyhzBa0AGQgVOGMFvSwafwFGrWijipXaFfkcccx4lTcwvhJG+8s5K1sV1B0n1pz8mDl2swh1ZUBoOMNnQAZChJU+FMC6ck1bhbYtbIu1uIVdnkJUtuwTdv70Z0zRt9G+/EqbnUhG5Skrzbpjmk3tABkIGTJLKW556N21uJ0tBYhgYq7eSs59CR2US9/SzCUE9ekqSd/ffA4u4kaQw7aJV2ja/w8U4trESigW8miQ5ARO4E8E5VfaRfg4vIZQD+AsAZAP6rqt7Qr7HIaOANGpC+cs2bNPUxcksET08tTRgPClvt+/LT2Hv4XYHp/MSE7wx/MrfowbO8TpuOY7hI2wF8EsDdInIrgJtU9SdlDiwiZwD4KwC/BuAxAN8QkS+p6nfKHIfUlyI7gVDbJxRh85U7afH7rNV3nnaJvZLlUMyQp638DZtvzPjHwka28k9a6bOBe7NIdACqepuI7AHwAQD7ReRTAObd9zf3OPZFAB5W1e8BgIh8DsAbAdABkCixBuQxYqWVtorOG0u3nUY/Vv9eeC3U7skjO5Ek4RxiuYNw19QLWSt/hpCGi6z/2n4C4ASAcQBnwzmAEnghgEfd58cAvDK8SEQ2A9gMAJOTkyUOT6qmiNHYsenihYSrXwHHDoB5zZuQrBV9aFx70ehPwj8vz6o9D6Hj8DuYE6cWy0z9SWZg6S4sKblORpO0HMBlAG4G8CUAr1DVH5c8dkyAdsn/aqq6E8BOAFi/fn0fNuSk7tjpWjOcoWhbSCzuHWu44pOtSSeE+xECyiKcY+gkpqcmlszXi7nZ770zCQ172XST2CfVk7YD2AbgTar6QJ/GfgzA+e7ziwA83qexSA3pxWiYBo9/VugorDomdioXwMK14UGpcJxBOYFl0qrZt/mETdaBpUqd/ndA567KdP5j8w//rmnAm0laDuBX+jz2NwBcICIvBvB9AG8G8Pt9HpMMIRa2MAO/bnJiicHasmvfkgbnWeGUPBr7MePpQ0x5evnmZV7R4azssJnXEgr7CpgjsJLOpBPFg9Lmp+MYLiprCamqcwDeBWAvgAcB3NbH3QapMTs2Xdy14fD3hnHwPKeA81xnWKLVnIbtHnrBnmmtGdMqlGLs2HRxx07It2cEOpPZ4WG3pOfRiDeHSg+CqeqdAO6scg5keLCdQIgZNet+ZYli3yQltkovunL3xtiLqZ04Ndd1qMiv9m1nEeoGhaWtRloS3XZLSX0OCAF4EpiMIGF+oAixpGuMMG8QGv8kAbjw+aHjiP0cqoj6fr5J+INyjOuTJOgAyNCS1gsA6GxIbqEWH/ePGXufdE1b2Set9v1q3eZlY5pjCnvm+rkW2ZWkJW7LNPZ0IKMLHQAZGmKGKEtGOdTJyUoMh/IPWQ1VwtV8GqH4m6/Jj0k0eGfgw0O+pLOIjhIhIXQApHbkbSoeNi4H0HGflzwwfFgG6DTiyyTZQdh1fkcQXmv9g/3K/ODMbMd15mDmNXll7d/Hxrbnx8JT/TLwPN07+tABkJ4YhFEIDdHG7XsWkqOhNpBnx6aLF64NiSl5pq3203YEaSWWJtIW6xucRKxBTLhDoBEmZUAHQGqDrebN2B6cmV0QJ8uT1PUGPa2/rqdo9U7sejPYvhOZVfWYo7Fy0fAQV8yQm8OwZw2ifj8GD4eNPnQApDDWZtCHJGItGcsyGGYAbTXvq2MALEmuFiGP8feSE142okhf3SyRt5jzA5bmFWiESZnQAZBKSKtcMUM/PTXR8X2SsbY2h+aYej2ZG6u/jxn5cJzws4/923O86NrBmdnMg1nzio7rqtwJkNGDDoDkJlylWlLUG+Akrfm8mv9eCuHQkVYIyKplrEmLl07uR3gkdDTesHsRuiIyEjGJao/X4w8lLQjpF3QAZCCEDVxiTsIbd2BR0TK8x7BVtD8YVYY2j6ltdtMLwPccjp0zyIqrh8bfnA1X4aQfiOrwLDXWr1+v+/fvr3oajSerTDOUJADQcfjKG7kwzBIzuuE9wNLTtWldwHrBnIFPzKaNkyQUZ60mPeHfYVjS6qEDIL0gIgdUdX34+8rE4EgzOHl6Llpbb+JnwKJxy6r0CTtpheJs9qzpqQmsHB/rSN5mHdJKwmSVTUzt8NHjC/NdJq33sPGmpyawe+sGrJucWFAstXmEImuhiJsPrx2cme0Yi8af9AuGgEhhsgySzwdk1eCfODW3ZMcQCqKZobTqG1vxhzsGyxmYIqbPPRRR7UzS8UmL3yclc8Pkb0y2IinuT/E20m+4AyADwVbgMaNmq137ft3kxELIxIyircSBReexe+uGjpW+5QzCVfPho8dz9do1/LXrJicWVvo+Uet3BhbmAhZP/xat3PFjDkq7nxDmAEhfCFe8ZqRNshlYNORJ8XFf5x9T27Sm57HvAOCuay4HUKyyJhzLf7YqJD93e7dYDsLeyxvytMqovBIYhBQlKQfAEBAZCGHTEgAdVT/eIIY193mlHMLr8/S/nZ6aSDx0FY7tdymhomeSkS8KjT8ZJHQApGfSDnWlGUOL4x86Mtuh2dNtwjYkrL03DX+/qvfOxodhQueS1JQlbdwYaZLNNPxk0NABkL6RVtZov++m1j4voayCEQsZmayD4Us/Y2GZtHBN0eQttXZIVdABkK7JIxcc1rn7TlpJ9fq2Su8H/qCWzc9W7BaSWjc5ERVxyyLMbdCwk7pDB0BKJ8kxAPmVN8uSQvBnAYBOATmgM/cQk5aOVQ8V1fNJgnr7pGroAEjX+LJI/zlmEC2k0gthhU7e/r0xrILISzPHxsuzmg8rgyjfQIYFOgBSiFgZY5KRDI1rVjvGGGl6/UWMv2+/aIRJYluJZ0lLh44vvD8v1NsnVUMHQAbCusmJqNaPGfgwNh+jrLCQX7EfnJldqNf3Oj5+vKx8BA05GVYqcQAi8iYAHwTwSwAuUlWe7qo5afHqmOEL4+ShcQUWD3MByavoPE3ZV46PYffWDR2tIo1wdxL+HlisSPIJ4bQ5JO18uoFOg1RJVVIQhwBcCeDeisZvDKbJUzVm/E26YeX42EKZpTkJo2gFkMkxeKmIvOGh8CCXOaJwDl7eISlnQOE2MmxUsgNQ1QcBQKRPtX4kF0VWn2lhjtj9tjOwVbo3zrZCj9XL5zH+tnOwuP7G7Xui+YXdWzcs2ankad8IpO88ijRoT7qGFUCkDtQ+ByAimwFsBoDJycmKZzM81MXAJK3E/arbtHqAuApniNXo2zutm5yIqn1u3L5nSYP1mCaQNXT3Gj5h1y9W9pBRpG8OQES+AmB15KttqvrFvM9R1Z0AdgItMbiSptdo8jiHWEVP2jVJO4MkeefQkM7r0vp7O5TlD495/Encjdv3LPl+XjvlmNPkHcKdgZePiDWDCd81JOvvmIljUgf65gBU9dJ+PZsUJ2/ooyzCRKkRE2pLUta0Od91zeULp4iT5BnCKiMf2w9zIJaQTurGVZawGyF1p1I5aBG5B2BMNKIAAAzYSURBVMB781YBUQ46P2lyzOE1MYMXKmSGksihIqaXdo6pZnr8NUCnA7D4flbYxj/Lv4NVAk1PTSTG3Q3vHNLUPXuBDoTUgVrJQYvIbwH4SwDnAtgjIver6oaM20gBwsNKg9an8avvUBAuXJX7lbuFfrwDiom3JTkEE3VLS9KmvXs3Xbho5MmwUlUV0BcAfKGKsckiadU7thJPS4KGOYBDR2Y7+tr6+L+dBbCSzbAZTFJzeN/43XICwGJsvohxj+0IktRKy4JOgdSZ2lcBkd4IJRB6MUimnVPkeu9E1q5eFa3WCROuaZ3ALHQUVgKlhbH6tUKvS6UVId1CB0AWCA1amGy98qa9HYnYsJLFjLLPCfg8wOGjxzuuATqNpe+ta5hWvw8nxU7x9vKeAPvwkmZCB9AQujFstlq3OnkzvBbisdLLsJmKr+mPrfhDwoohf1DMVvw+XJOVsB1UiSVLOcmwQwdAFoipePq2jSEm4hY2fPFhHMP0evIYS3MgfoXuq4Z8JVI3RpeGm5AWdAAkmsz1XbvS5BnmtVMW2ap4/PfAYrzeh2ySDpsZPnQUNmvJY/wHZdjpQMiwQgdAlhCu4P3hq7SafF96mkcCwojV59vvwwNsPn/g7+1lJ5AGdwlklKEDGFGKCJWFVSwWqvFhFm+kfZLXSzsknf61e4wrb9q7JJZv44bESkNZq09IOdABkAWyKmr8d2a88zZpOXl6LrEtZNIhtbBCKKwI8gJuZRt4lniSJkAHUHOKGp60huzhc9J0cPx3oVKm/TN2j5eA9lhuoNsevknNWrLo1pAPWjuJkCqgAxhxslb1Pulrp3ezEqxhS8U0o+oPc/ndQqwBTOz+WM3+ll37+l7JE/6dceVPRhE6gJrS7co1KYGa9pxY2Wb4PKv0iQmoJY1vhPX7QLqej1GGVENRR5EkokfIKML/ukeAWDllrHTSfg7LML3x9/H40FiaFEQRoxpKNIQni4F48tfI6r7Vzcq8yK6hnzpBhFQNHUBNKbpytdO6dl1sxR1KOOQhdBLhOFnzD8kzfmz3U9ZOIO91TPqSJkAHMMSEhjI0zvZ9Vsw+LVHs77dnHDqymCsoQvjctJV/SNr5gV4T5DTypKnQAdScItUqJ07NdTgBb/h7GT8MKYWlnFmtIYtQ5rN6gU6BNAE6gCEmZpz9qtx+DhOasTi63wWEIR4vEeETxr5pS6y+354VG7MXykiQlz0nQoYROoAhJ+y6lVTnH1LUCPqkbYhJPHsJiCJ9A/z1DM8QMjjoACqiTANnxjnr8FLeOvuk+fnmMjFdIC8CF4afyn7fXp5Jp0JICzqAESEtIZvWLatInXu4qg9bRsbmEDsrkDZHrvwJGRx0AAOm7FBHnueljWESC16+ISlXYKTV9idVGPXDoNNJENIbdAAkkywnE9MQIoTUHzqAAVP2yjjteXkPVGWdts0zfjhmLw4h72EzQkhv0AHUkF6cQ3hvmBjuxth347R6PaRFJ0BI/6nEAYjIDgC/AeA0gMMA/p2qHqtiLlVR1DAWrXG33yWViPYyp7R5+rMBRUk71EYIKR9RzdnRo8xBRV4P4KuqOiciNwKAqr4v677169fr/v37+z6/GIOoTgmNaJHm57HqHl+maVU6/Vi9x8Y3rH9A3nH9obYi9xJCkhGRA6q6Pvx9JTsAVb3bfbwPwO9UMY8sqi5JNGOYpbmfh36rWtq8Nm7fA2BR6rlIY5WiOxZCSG/UIQfwdgCfT/pSRDYD2AwAk5OTg5rTAhaGGMQJ1VDOOVT07PberLJOo4wSVZOE6LZ/byhLQQjpH31zACLyFQCrI19tU9Uvtq/ZBmAOwKeTnqOqOwHsBFohoD5MdQmxhGQvz+nGWWQlawe5Oykylp0W9g1kisKVPyGDoW8OQFUvTfteRN4G4AoAr9MqEhEFsMNSgwpL9DJG3gNgafdmyUXngY1UCKk/VSWBLwNwM4DXqurTee8bdBK4W2niWDI3771lPrvbXUJo8PvxHoSQwVGrJDCAjwIYB/BlEQGA+1T1qormkouqjF4voZ5u59zLyp8QMjxUsgPolirLQLuh28NQecI4dc0BEELqR912ACQD6uMTQvoNdwA1IO0AGOPvhJBe4Q5gyOhFoI0QQvJAB1AD2AyFEFIFdAA1h86AENIv6ABqQJUrf+46CGkuy6qeAOmdLbv2sWafEFIY7gAqpMpST5aZEkLoAIaYLCNOo04ISYMOoELKNtSxfr+DGpsQMnzQAQwxsR4A/uwAwzuEkDToAAZMzBj32os37BzWzU6AENI86ABKpFv56F4JO4fxFDEhJA90AAOiX2EZxvIJId1CB1ACoXG/8qa9Cz1x0+rzrd9wv4w2nQEhJA06gAERW6mXeXiLxp4QUhQ6gBJIM+4xw2wrf1bpEEKqhFIQA2bHpovZMJ0QUgvYEKZCuPInhAyCpIYw3AEQQkhDYQ6gQrjyJ4RUCXcAIwjloQkheaADIISQhlJJCEhErgPwRgDzAJ4C8Aeq+ngVcxklKAJHCClCVTuAHao6raovA3AHgA/0e0CGRQghpJNKdgCqetx9XAlgeGpRawx1gQghRaisCkhErgewCcA/A/g3KddtBrAZACYnJwuPw7AIIYTE6ZsDEJGvAFgd+Wqbqn5RVbcB2CYifwbgXQCujT1HVXcC2Am0DoL1a76jBJ0bISQPlZ8EFpEpAHtUdV3Wtb2cBObKnxDSVGp1ElhELnAf3wDgoSrmQQghTaaqHMANInIhWmWgMwCu6veAXPkTQkgnVVUB/XYV4xJCCFmEJ4EJIaSh0AEQQkhDoQMghJCGQgdACCENhQ6AEEIaCh0AIYQ0lMpPAhdBRJ5G69xA2ZwD4Ad9eO6gGYX3GIV3AEbjPUbhHYDReI9e32FKVc8NfzlUDqBfiMj+2DHpYWMU3mMU3gEYjfcYhXcARuM9+vUODAERQkhDoQMghJCGQgfQYmfVEyiJUXiPUXgHYDTeYxTeARiN9+jLOzAHQAghDYU7AEIIaSh0AIQQ0lDoANqIyHUiclBE7heRu0XkvKrnVBQR2SEiD7Xf4wsi8ryq59QNIvImEXlAROZFZKjK90TkMhH5rog8LCL/ser5dIOI3CIiT4nIoarn0i0icr6I/E8RebD939KfVD2nbhCR54jIP4rIt9rv8aFSn88cQAsRWaWqx9s//wcAL1XVvjeqKRMReT2Ar6rqnIjcCACq+r6Kp1UYEfkltJoFfRzAe1W1uz6gA0ZEzgDwfwD8GoDHAHwDwO+p6ncqnVhBROQ1AJ4FsCtPq9Y6IiIvAPACVf2miJwN4ACA3xzCfxcCYKWqPisiZwL4GoA/UdX7yng+dwBtzPi3WQlg6Dyjqt6tqnPtj/cBeFGV8+kWVX1QVb9b9Ty64CIAD6vq91T1NIDPAXhjxXMqjKreC2C26nn0gqo+oarfbP/8DIAHAbyw2lkVR1s82/54ZvtPabaJDsAhIteLyKMA3gLgA1XPp0feDuCuqifRMF4I4FH3+TEModEZNURkDYCXA/h6tTPpDhE5Q0TuB/AUgC+ramnv0SgHICJfEZFDkT9vBABV3aaq5wP4NIB3VTvbOFnv0L5mG4A5tN6jluR5jyFEIr8bup3kKCEiZwG4HcC7g13+0KCqP1XVl6G1o79IREoLy1XVFL4SVPXSnJd+BsAeANf2cTpdkfUOIvI2AFcAeJ3WOMFT4N/FMPEYgPPd5xcBeLyiuTSedsz8dgCfVtXdVc+nV1T1mIjcA+AyAKUk6Bu1A0hDRC5wH98A4KGq5tItInIZgPcBeIOq/rjq+TSQbwC4QEReLCLLAbwZwJcqnlMjaSdPPwHgQVW9uer5dIuInGvVfCKyAsClKNE2sQqojYjcDuBCtKpPZgBcparfr3ZWxRCRhwGMA/hh+1f3DVslEwCIyG8B+EsA5wI4BuB+Vd1Q7azyISK/DuAjAM4AcIuqXl/xlAojIp8FcAlaEsRPArhWVT9R6aQKIiKvBvC/AHwbrf+nAeBqVb2zulkVR0SmAdyK1n9PywDcpqofLu35dACEENJMGAIihJCGQgdACCENhQ6AEEIaCh0AIYQ0FDoAQghpKHQAhHRJW3Hyn0Rkov35+e3PU1XPjZA80AEQ0iWq+iiAjwG4of2rGwDsVNWZ6mZFSH54DoCQHmjLDRwAcAuAdwB4eVsJlJDa0ygtIELKRlV/IiJbAPwdgNfT+JNhgiEgQnpnI4AnAAxl8xTSXOgACOkBEXkZWh3AXgXgT9udqAgZCugACOmStuLkx9DSmj8CYAeA/1TtrAjJDx0AId3zDgBHVPXL7c//GcAvishrK5wTIblhFRAhhDQU7gAIIaSh0AEQQkhDoQMghJCGQgdACCENhQ6AEEIaCh0AIYQ0FDoAQghpKP8fIPVHspjC7+4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[0,:,0],Y[0,:,0],label=\"data\",marker=\"+\",color=\"steelblue\")\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.title('Plot of data samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the MINEE model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from model.minee import MINEE "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_size = 100        # batch size of data sample\n",
    "ref_batch_factor = 10  # batch size expansion factor for reference sample\n",
    "lr = 5e-5               # learning rate\n",
    "\n",
    "minee_list = []\n",
    "for i in range(rep):\n",
    "    minee_list.append(MINEE(torch.Tensor(X[i]),torch.Tensor(Y[i]),\n",
    "                            batch_size=batch_size,ref_batch_factor=ref_batch_factor,lr=lr))\n",
    "dXY_list = np.zeros((rep,0))\n",
    "dX_list = np.zeros((rep,0))\n",
    "dY_list = np.zeros((rep,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "load_available = True # set to False to prevent loading previous results\n",
    "if load_available and os.path.exists(chkpt_name):\n",
    "    checkpoint = torch.load(\n",
    "        chkpt_name, map_location='cuda' if torch.cuda.is_available() else 'cpu')\n",
    "    dXY_list = checkpoint['dXY_list']\n",
    "    dX_list = checkpoint['dX_list']\n",
    "    dY_list = checkpoint['dY_list']\n",
    "    minee_state_list = checkpoint['minee_state_list']\n",
    "    for i in range(rep):\n",
    "        minee_list[i].load_state_dict(minee_state_list[i])\n",
    "    print('Previous results loaded.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuously train the model. The following can be executed repeatedly and after loading previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEICAYAAACtXxSQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhU5dn48e89S3YIS1gCaHHBrVbRl1prbat1Lb9au7xt9bV20com4I5IIBmSoCIqisimYtVqrXvtolVbl9pWFK1WRauoKEgSAoTsyyzP749zksxMZpKZZCaTTO7PdeXKzFnvnJnc88yzHTHGoJRSanBzpDoApZRSfafJXCml0oAmc6WUSgOazJVSKg1oMldKqTSgyVwppdKAJvM0IiIviMivUhzDoSLybxGpF5H5MWzvEZHf2I/3F5EGEXEmP9L0JSLnicgzqY5D9S9N5oOMiGwTkWY76VWJyN0ikhfnMSaLiBERVxJCXAC8YIwZZoxZFc+OxpjPjDF5xhh/EuJKS5FeS2PM/caY05N0vpQXGFRkmswHp7OMMXnAscCXgcUpjifYF4B3Ux1EsCR9aCk1oGgyH8SMMZ8DTwFHhq8TEYeILBaRT0Vkl4jcKyL59uqX7N/77BL+V0XkYBF5UURqRWS3iPwu2nlF5Lsi8q6I7LNLaofby/8GnAysto97SIR9D7DPUy8izwIFQes6Spkico6IbA7b9zIRedJ+nCkiN4rIZ/Y3lHUikm2vO0lEdojI1SJSCdxtL18gIhUislNEfmWf6+A4jneFfS0rROSXQXFli8hN9rWuFZGXg/Y9XkT+aV+rt0TkpG6u6wQReVREqkXkk+BqKhE5TkQ2i0idHd/N3byWvxCRl4P2NSIyR0Q+tK97mYgcJCL/so/3kIhk2NuOFJE/2jHU2I8n2euWAV8Pen1X28sPE5FnRWSviPxXRH4cdO7pIrLFPu/nInJltL9f9ZExRn8G0Q+wDTjVfrwfVim4zH7+AvAr+/EFwFbgQCAPeAy4z143GTCAK+i4vwWKsD7gs4ATo5z/EKAROA1wY1WrbAUywmOIsv+/gJuBTOAbQD3wm/C4gBx73ZSgfV8DzrEf3wI8CYwChgF/AK6z150E+IDl9nmygTOBSuCL9rHvs891cBzHK7X/5ulAEzDSXn+7/XdPBJzACfZ5JwJ77O0d9jXbA4yJcF0cwOtAMZBhv24fA2cEXbfz7cd5wPHdvJa/AF4Oem7sv224/fe3An+1z5EPbAF+bm87GvihfY2GAQ8DTwQdK+T1BXKB7cAv7dftWGA38EV7fQXwdfvxSODYVP8PpetPygPQnzhfMCuZNwD7gE+BNUC2va7jH83+Z50TtN+hgNf+h4uUAO4FNgCTejj/EuChoOcO4HPgpPAYIuy7v50Uc4OWPUCEZG4//w1QbD+egpXccwDB+kA5KOg4XwU+sR+fBLQBWUHrN2InZ/v5wfa5Do7xeM1h12sXcLz99zcDR0f4e6/G/gANWvaX9sQZtvwrwGdhy64B7rYfvwQsBQrCton0Wv6Crsn8a0HPXweuDnp+E3BLlNdsKlAT9Dzk9QV+Avw9bJ/1QIn9+DNgJjA81f876f6j1SyD0/eMMSOMMV8wxswxxjRH2GYCVrJv9ylWIh8X5ZgLsJLaq3YVygVRtgs5rjEmgFUymxhD3BOwEkNjWFzRPACcaz/+P6wSYhMwBiupv25XX+wDnraXt6s2xrSEnXt70PPgx7Ecb48xxhf0vAmrhFyA9U3mowjxfwH4Ufsx7eOeCBRG2XZC2LaL6Hy9LsT6VvS+iLwmIt+JcIzuVAU9bo7wPA9ARHJEZL1dZVSH9SEyQqL3MPoC8JWwuM8Dxtvrf4j1zeRTu3rtq3HGrWKkDUPpayfWP1q79lJxFRESrzGmErgIQEROBJ4TkZeMMVsjHPdL7U9ERLCqez6PIaYKYKSI5AYl9P2xSo6RPAMUiMhUrKR+mb18N1YC+qKx2g0iCT9mBTAp6Pl+QY9jOV40u4EW4CDgrbB127FK5hfFcJztWN8EpkRaaYz5EDhXRBzAD4BHRGQ00a9db12B9S3uK8aYSvva/xvrg54I59sOvGiMOS1K3K8BZ4uIG5gLPETotVcJoiXz9PVb4DKxGhzzgGuB39mly2oggFVnCoCI/Ki9oQuowfqnjdRF8CHg/4nIKfY/6BVYdbD/7CkgY8ynwGZgqYhk2B8aZ3WzvQ94BFiBVZf9rL08ANwBrBSRsXb8E0XkjG5O/xDwSxE5XERysOqm28/Tm+MF77sRuNluwHTajZCZWNVEZ4nIGfbyLLsxdVKEQ70K1InVaJttb3+kiHzZjuenIjLGPt8+ex8/EV7LPhqG9cG2T0RGASVh66vCzvVH4BAROV9E3PbPl+3rnCFWn/d8Y4wXqCPye0olgCbz9LURq5HvJeATrNLjPAC7qmIZ8A/7q/HxWF0cN4lIA1Zj2SXGmE/CD2qM+S/wU+A2rFLpWVhdJdtijOv/sOqH92Ilint72P4B4FTg4bBqjquxGl5fsasDnsMqUUZkjHkKWAU8b+/3L3tVa2+OF+ZK4G2sBtq9WA2vDmPMduBsrOqSaqxS7FVE+L8zVt/6s7DqqD/BurZ3YjVQgtWA+679+tyK1RDcEuW17ItbsBqMdwOvYFU3BbsV+F+7p8sqY0w9cDpwDta3tko6G54Bzge22dd0FtZ7RyWBGKM3p1BDj1jdKd8BMsM+JJQalLRkroYMEfm+/dV/JFbp8Q+ayFW60GSuhpKZWNUdH2HV3c5ObThKJY5WsyilVBqIuWQuIvuJyPMi8p7dD/kSe7nHHqb7pv0zPXnhKqWUiiTmkrmIFAKFxpg3RGQY1iiy7wE/BhqMMTfGetKCggIzefLkXoSrlFJD1+uvv77bGDMm0rqYBw0ZYyqwBl5gjKkXkfeIbdRfF5MnT2bz5s09b6iUUqqDiEQdMd2rBlARmQwcA2yyF80Vkf+IyEa7p0CkfWaINevb5urq6t6cVimlVBRxJ3N7NOGjwKXGmDpgLdZQ5qlYJfebIu1njNlgjJlmjJk2ZkzEbwlKKaV6Ka5kbg/ffhS43xjzGIAxpsoY4w8aEn1c4sNUSinVnXh6swhwF/CeMebmoOXBM8B9H2tUnVJKqX4Uz6yJX8OaZ+FtEXnTXrYIaya3qVgTM23DGpihlFKqH8XTm+VlOqfBDPbnxIWjlFKqN3Q4v1JKpQFN5kopFSef18tDy4upqdre88b9RJO5UkrF6S933MCWZgcPrr4z1aF00GSulFJxamn2AtBqot0atf9pMldKqTSgyVwppdKAJnOllEoDmsyVUqrXIg29SQ1N5kopFadIKfz9F1/g1uLltDY08uLTrxIIBKiu2MuWNz/ql5jiGc6vlFJp641Nz/DGX57hV8Ux32cnxOPP/ZNWZxsbb7+fquadtDa38sq/XyLg8OKZ6klssBFoyVwppYC/PbmJHYE8Hr0n4izeMWtsbQagtraOgMObiNBioslcKaUAr1jp0NvcnJDjWRPN9h9N5koplQSfbt/Wr+fTZK6UUiFiu8l9T+p9uxJynFhpMldKqYRKzIdBvDSZK6VUktXV1bF79+6knkOTuVJKBWmsruV3G5aFLLunfAG/Wb4gpv1NhIL5zTffzOrVqxMRXlSazJVSQ8bdNy3iXy88HXFdwP693Z/Hezu93FhUzrZt7wHwiS+Hrc05XfYxGGoq+mdQUE80mSulhoxP6zN46bnXI67zukL7hDe4fTz/8N0hy174/Qbrgd3rsNZpuHX9fWx64u6OmvIWf2MiQ46ZjgBVSg0pza7eD+RprN0bukCsFP70GzswTj8A/n4cKBRMS+ZKqSHv7msvj2v7tjZ/yHPj8EfZsv9oMldKDXmftg2PuLxuX9d6coD3WwbOHYbaaTJXSqkoatxdl21/b3P/BxIDTeZKqQGrtqaGzz7bmpBjNTUmpmHykfufSshxEk0bQJVSA9Zvb76eSmc2047an+/84IKUxODxeEKee40DSH0deTgtmSulBqzdWHXZOz+y+nuv88zl7c0vJez4zU1N3FUeX+PnQKXJXCk1KKxbOpdKCnjxiWd7tf8LTz3SZdmqa1ey3Re58XOwiTmZi8h+IvK8iLwnIu+KyCX28lEi8qyIfGj/Hpm8cJVSQ5XdpZsWk8Fd113O+tJ5XbZ5YE0pDXV1Efd/+80dIc/vKL20T33Oe+ONN95I2rHjKZn7gCuMMYcDxwMXi8gRwELgr8aYKcBf7edKKZU021uHUxEYHbLszuVX8MGuAL9ZURZxnzYJnTRlb9uYpMUXzZNPPpm0Y8eczI0xFcaYN+zH9cB7wETgbOAee7N7gO8lOkil1FAV+3SygWarlN0YGNjVJrW1tUk5bq/qzEVkMnAMsAkYZ4ypACvhA2Oj7DNDRDaLyObq6ureRauUGlI6UrkIFRT06Vh+p6/P8STCa6+9lpTjxt01UUTygEeBS40xdbHe584YswHYADBt2rTUzN6ulBqUHBK93Hn3DQvw+zvrvuvdfjweD1PG5vBhVQuFZh/Hnvb/QvdZsZBmV1bccWz/sB6/RB4VmmpxlcxFxI2VyO83xjxmL64SkUJ7fSHQv/dKUkoNWdu3fcCnTTnsaM3vsq6+YjdIgArHcP7+9D9C1n3aGH8iB6h0umnt50bTWMXTm0WAu4D3jDE3B616Evi5/fjnwO8TF55SSoXrrA3488Z1Me3R5Ej/XtjxVLN8DTgfeFtE3rSXLQKuBx4SkQuBz4AfJTZEpdRQF7U613QmadNNY+lQqNeNOZkbY14m+CMx1CmJCUcpNdSt81xMwJfPLlcGTjtFba/LDNoicmpu8haA20TcbqA0fiaTzs2ilEqYNeWX4mx1MbPsxt7tv3gRu1xjus1MDe7I86LUukOTfJ3JxxoeMzSkf0WSUqrf7GsZQ4Uzj3XXRZ7vZH3JPDweD3ffWgpA5eefscYzh5o9e1i3ZAG7XBlxna/CkRd1XZNr6CRy0GSulEogvz3K0gQiV4XU+McD0FJTBcAT61azi7E8dGs5lc5edPmLrWf0kKDJXCmVeBEaLB+9dw0tdre+fb4JrC++jNaAlcD9gex+DS+VfL7kfGPQZK6U6hfb/9s58rvV5aXC0bVv+FCwffv2pBxXk7lSKuEi9SRsI8J9M2McQa56pslcKdXFA2uX8/i9sQ3IaVe2pDzhXQCHQpfCRNGuiUqpLj6oasbp9/L9OPYJTrymm7lUgtW47IZSGQrDeixeb3KmA9BkrpSKqLtS8T23L6WpooXjv/dtPnjzVWp2VQDDOvdtaY7rXN5AJpFqYVTsNJkrpeKybskCqxuhK5PXnnyEKjMWv3NYyDZ7A4Udj9cXX06FYzi4ox+zSdzAwJzAarDQOnOlVFyC+4MbHBFL8AFH57IKR883ixioMxEOJprMlRriVnnmsq5kbo/brfHM4fZll4Usa/VH715454qiPsemYqfVLEoNcfsC40NK0uF+s+56avfuppqxOFpcEFQS3+uO3rVwR6Ob9SXzQUYlNF4VmZbMlRri2hP5zu3b+fiDd7qs/+RzH9Vt1hwogTi7ClZoIu/CmOT03NGSuVJDyIabSjD7Gpl+4Vx+d8d95EoDuKyh9E+sX8suVwannrAtZB/t6z04aDJXaghp3ivUuPP4090raHCPoYHOOVHaZyz85IP/pCo81QdazaLUEOJt78wdiD5DYU1Fa8fjNZ7ZyQ5JJYgmc6WGoEpnbtR1wY2auxjXH+GoBNBkrtQgdmuRh6XFpTFvr9NapS+tM1dqEKtxAwRi3n7ozIAy9GjJXKk0sW7xQjwlnlSHoVJEk7lSaWB9yTwqXVk91qNoNUvqJaufuSZzpdJAvW9sTNv5NZ2nLU3mSg0RKxaX0+TWAUDpShtAlUqBVcWzcYhh7tLu7+azvmQeJpDFrLIVXdat8cwBrBJ5g9vfsXx1yWy++/Or+PM9N1MpBRQG6jCOFhpdsZXe1eCkyVypFNjriK3/doWMjnrThl1ETs67ZRx/2LiRaleBdQzHcKDnaWjV4BZzMheRjcB3gF3GmCPtZR7gIqD9ttuLjDF/TnSQSg01a5dcQ43J6/Ifurr8CnzeJi5durbb/atdWk4bqAZCA+ivgTMjLF9pjJlq/2giVyoBqpyZtIXdsGF98WXs9g1jn+ioTNVVzMncGPMSsDeJsSiloli/5AoqHJ03gtA5U1S4RPRmmSsi/xGRjSIyMtpGIjJDRDaLyObq6upomymlIqgIu8emzpmiwvU1ma8FDgKmAhXATdE2NMZsMMZMM8ZMGzNmTB9Pq9Tg9dKzT3c8fnPzyx2P1y5eyPLFy/jo/fdTEZYa5PqUzI0xVcYYvzEmANwBHJeYsJRKPzeXzOXZJx5k8wuvdix74o/PsbpkFgBVriyaXV7+8MCqVIWoBrE+JXMRKQx6+n2g6z2nlBpCbi3ysGbxoi7Lb1kyizop4B9vvk+dO3RiLLfJZmVR58yH+xzaH1zFL56uib8FTgIKRGQHUAKcJCJTsSZj2wbMTEKMSg0a1iyGGR3P1y25mkpnNoWSw75u9qt1xz7zoRrcUn4PUGPMuREW35XAWJQatFZ75uAPBMAxvmPZyiWzqWci4LMH7iiVPDqyQKk+Wu2Zw27GhlRarvHModY5Duh5LpTgLocq/TmdUYb09pFOtKVUH6xcMtNK5GEavRNSEI0aDA444ICkHFdL5kr1wpqS2TiNm1pnYcT1jTo7oepnmszVkHZ7yWycOJm1dHXM+6xcMtuqQtGpwdUAoslcDWnVUeY5WV8yn8ZAG5eXrQtbPo8ccqjtj+BUWhoIE20plVZ+e1dnafyOW8s7Hq9dvpgKGUWz2Q+Px8OqxSUd6ypkdJeh9UoNBFoyV0PW1m37OuYK/7zGh8fjYQS1OFxOIA+f07rhw16X1qeoxBFJzvtJk7lKWyuK5+ALTCTHsZ1LStdxS8l8sg3MLLWGy/udXRspswMBWrTtUiVRygcNKTXY5JPJTpeXHK89kMc3mgq3YW3JxVRJ5MneKhwjGe9vBCcYCR2Vub74MtA+4WqA0jpzlbaqA9Zt0wI4WLdkAS1ilV2iJfJ2lc7cLsvWF1+ug3vUgKbJXKWFt17ZxMrFszqeryuZi9e+U0+tO0ClM4eWsDv3xEOH46uBTpO5Sgt/+8NfqHWNx+PxAFApBakNSKkoCgqS897UOnM1KN15Syk79gUoDNQSEC+17s5/kBVL5oBTp5FVA5Pb7U7KcbVkrgasVcWzO27cEG7f7p0AVMgIHCYjZF1WQEvlauDS3ixqyNnrsEZn3lrkocYNLr+bxWVFAHR01RWDCYS+jU2S+vEqNZBpyVyl3KqSWawvvjxk2V+febzjcY39rdTntBowr1+8jHpn57zhAROazHWQjxqKNJmrlNtnJlLhGM5zTz3Wsezv/3wr6vbhvVL2SteuhEoNNZrMVcoFHNaw+RG5Vv/vmz2XRN329pLZXZa1l9iVGgySNZxfk7kaMN547mEA6hgZdZtosxwqNdRpA6hKidUls6gJTGKcVIPDSt47HaPsUnn0ZK6UikyTuepXNywup8nlI983Ab/bx86wxG28BeD2pyg6pQYvTeYqadYvuRK/NDKndC1LS0qtiavsd1ytOxBxn3pN5Er1itaZq4S6sWQ+z/zhQQAqnHnssvuKh89AqJRKLE3mKmGWF82kQUbxz9ffD+k3/uyfH05hVEoNLFu2bEnKcTWZq4RYVTKTZnfnneqDZxn81yvvpyIkpQakxsbGpBxX68xV3DauW0H9zo+ocYxnnNnF7KVryAjkd9yCLVx7P3KlVPJoMldx27GzhYDDGk6/OzCB8iXL8DlzUhyVUkNbzNUsIrJRRHaJyDtBy0aJyLMi8qH9WzsIp6FVJTNZFTR7YXBJ2+/06QhMpQaAeErmvwZWA/cGLVsI/NUYc72ILLSfX5248FQq3VYyG3dgGHudVl34bUXFZLhqQEanODKlVLiYk7kx5iURmRy2+GzgJPvxPcALaDIflP7xzHO89uI/2Oc2OAJOikuXsEfGhdSD73E7AE3kSg1Efe3NMs4YUwFg/456excRmSEim0Vkc3V1dR9PqxLt1b8/wj63NWl+wOHvuP2aUiqxBv3NKYwxG4ANANOmTUvOX6PicnvJbJyBPIxxUusa3/MOSqkBq6/JvEpECo0xFSJSCOxKRFAqOV5++hm2/PMpwM/uwDjaXOOididUSiXHUUcdlZTj9jWZPwn8HLje/v37PkekEm550RyGuQ313kk0u/OthQ7tgaJUKgwfPrznjXoh5mQuIr/FauwsEJEdQAlWEn9IRC4EPgN+lIwgVd9kMo5dGHBrAlcqXcXTm+XcKKtOSVAsKkFWLJ7BcEc2FY6RjPc3sc+tA3qUSnc6AjQNNbom0D77Q6WOzFRqSNCJttLArcXW6MybSi7hxqLyFEejlOpOsu4BqiXzQW598aXUOMazanEJjTKGgNuX6pCUUimgyXwQWr54DsOcIP7hVLlGALDXJYAmcqWGKk3mg8zyootodk+kGfTVU2oQ0mqWIe7hBzdS8d4mmt0TUx2KUmoA0mQ+SLz7/mcghT1vqJQakjSZD2A3FZVR7/YzztcKrsxUh6OUSgCtZhlCVhbNpNZdCG7reZUmcqVUD7Sf+QCzsqjUSuRKKRUHLZkPAMsWXYQ3w27YdKc2FqXU4KQl8wGgI5EPQoX++j7tP87XkqBIbKb7+shsn35aqvSkJfMUWlk0k3rHfgNiTvEMn5s2V8+zKuZ6XTQGjzLtRVtOAVXsZly32zgCzo4bR4831VTKmJhidPtdeLvZZhiNNJMRf9BKDXBaMu9nK5fMZEPxfJYtXkatu5CAs39GbRZEuW9IQaAKj8cT8Y2Q73WQEVaSvWrZ4pDn9RKaODNjKPnO9azteBzts6C4dEnH41lLbyffX0leTvffAkZ7A2TL9h7PH4nT37/lmjzvAPgEV2lFk3k/q3UWstMxqtvSY7yGBSWGCYG9ZLVVdDx3BFzkeyuY61nTsawwUAfAGJ+PuaVWYs1mR5fjZkkjJ3/zeHK9VRHP6/C7uHLpKgi6CWC+NPTpb4nmsrJ1zF90fbfbZLhquLxsXbfbGGdNIsPqtWxJcPWSGjSS1TVRk3k/WeeZ0+ubJI/zNXe7vpXteDwe8v2VzChdxcJr1wMw3OuguHQxly1bH7K9X5oZZ6rxZnZ+oFxSvpbstp04Ap0lVHE28dXTTuGqZWsZ62vrWD4hUAvAWEdlyHELAzU4R2V1G2t4SR9nvb1vbcg5kqUtEPn2s8n591Kq/2gy7yeVjI1vB7shL8/rZHb58o7Sd65vpzWIKMiiZRsAq/TaLqutgpwo1dIGmL30di5dsjxk+dXXbiCbKsaYXYwKVDFr6W0d6+aUX9vxYTSjdCUej4dZS1eH7D/qsCOZeVlJt39WtvGHPJ+19HZyvJXMLF3JnPJru923Rz00fgKc+t0Lo6zpXTof7fP3vJFS/UAbQJPoucee5PXNf6C5F71V8l2tTD36RE7+7nQAskf5ad2zk6uu3dCxzW0ls6O2nbaXziPKLoi66qrS26Kui9U4U40jkE2FMy9k+Xh/I/VilcQzfW6G27fQWLCs+6qRvsjzOmlwdybcL335y/z994+zyxVbI+h4fyOVztyo6zMcTcCwvoZpn6uZSmd2Qo6lhh5N5kly1+3XsbPSjz/GRO70uxhDHc0ZLi5b0rWEOucKT5dl85au7bIsFhcvLOvVftEUBuqpcAxj/3EHA1apH2DdkqtDktOsshUdj68pL0poDACIVYUy0ifUuKzHw5w1+HyjaQlqo3BKK4T1aMkKGBrsT8axvraOZH/QV4+n8tW344vDCOP8TVS5IifmMb42xFkDYT16vI5aoPfJ3BFwgSFqo3qsPZbCFQb2UuEY1eu4Qo6V1URFy9C++5UO5x9EVhTNp9E9qtsuh06/C7/Th9PvYknZ4ugbJtDYQBVOuq/T7o2ZZTdFXD6rbDkrimbS2F8jWu1qlkxpAKKXpiPJc9Rzpaez2qm37RvtTjnnfB545JGI6y4uv5bbSmaHLMvyuZlXvq5P5y0uXcxtRcXscUauPY2nTrXA58cYYY/bgZFAr2MKN3PhDX2+tioyTeYJtm7JVVYij6LQ38DMshv7MaJOc0p7V5Lvi6uWrWd98Xy8JK73TrBhXif17vjqrQOOZmCY1XggVuI66rRv9S6AiO2phkOOPJLMB39Pa5SS8Jk/mcf9Dz3U8bw/GmAzjZ94+tB855fn8OF7b/Hx5t1Ji0kljjaAJtD6JVd0W79aEKhKWSJPpZmlqzq6QCbatJOPYrg3vrfx7KVrGOavYoTPSqFuRx0nfPOMiNueNv2HURJ230w54ojEHzSK9r7/bonvA/WAg4/g9LPOTXg8+XG+XukmOzs57SJD+6om0G0ls6lwRm4Ik4ATj8eTtISWztw9DEL65qlnk+Haicvvps7EnqyuKIv9tRjtC61miGVgFCTlMyCGc3Yt44+S/u1bHz4AK/x6XbasuD/DGXAmTZqUlONqNUsCrC6ezR5HhH6ABjxLPf0eT7oY4d+F1/ScEucujd4bppXu++jHpNs6kN5XkBQGagAHFY78iOu7TJ3QZ2HX0khHo3GijfO1UOVKTPtMe/uS6p6WzPvo2msuYneERF4YqCfTuzMFEaWPS8vWcNWyyCVol7/n0vEYn4/53XwbcrmrkICDBtPDYCU7342MVPCX2OrrJdD1X21m6a0EpDXC1jA2yvQLvRGptN4TZx9vDh5wNXY8dsXwgdwdVwzjB5SWzPtk5ZJZtGWGdj3M8zq5ctmSKHuoRMl1VhEIRE6Enawkks0O8E0ik9B5W4KnOOit/aYeS8V/PupxO0HiqnaZ41nDiqLykGWJL6lHVhioY8oJJ/S4XbbPTXMMXR1HjWqlsb73c9GkorpqMEpIMheRbUA94Ad8xphpiTjuQHb/r9dS6xwfsizT5+bKZUnoPz2EFJo9tPq779Oc43Vx2bJVMR/zkvI+tlVEKBgGd6/b/ub8hPXD7s4wVxWNjO52mwzXHiDyoLDgPyPL58ZtAhF7As0svbkPUbbrbGe48Ioybigqpynog2hcZiOOlpSvjH8AABgOSURBVFYqJLHXrZA9VPRwjVLN4UhOhUgij3qyMWbqUEjkAB9uC518Kt/rSM5AmCFm5tLbmF++NOr6fH8FmVn7ejhKgr+W20XDDHtyrFGyN2T1zNJVjDPVvYopUaXO9ukeZi1dzbGHTiTX2305bWF5Ea4E9h/vyWFHTWKMo/O6zb5mBYFuqqhG+QwjvPG/jpl5sY8vOGhkgKMPGBn3OQYqrWbphbWLrwm5wfIon+HwY45KYURDx2Vl3UxTkGQO8eLxdD9zY3TJqSxoH9WZ5+xMlN899yJuLioN2qo9KYbGEF6VPdbfyi5ncu43+92f/AL4RUzbjvTC/GVLWV98GfuI3DjcV7leF+dfYg3WeytNBjElKpkb4BkRMcB6Y8yG8A1EZAYwA2D//fdP0Gn73xrPbHa5Ohs8x/pbmVN+XQojUuF8zgYgB4ezpxJ8jGIoIKaiiS7T5yaXNvb2cHaxS8AO/HR3X8I5ZZHfx16nBNeadOq2YbP3X/oz3KGNv3F9DMb4Qgxz1MZz1EEhUcn8a8aYnSIyFnhWRN43xrwUvIGd4DcATJs2bdC2aewKmk/D6Xcxp8yTumBURIlo2AwZch7vu9VOKMOd29nn3w+fs+dGwkiniO8DIvrWXzrpNMzzf2XiFw+l6r9d563vybzi23o9BL/bXkf2CNx2h47PJn9kAdN/0rtzxSMd+8ckpM7cGLPT/r0LeBw4LhHHHWjC39D9NaeKGlzaB83MXbqO4WIlT6fdNdHR5Y5GyS/XnPCtM5lVtoKzzv1VQo+b6Pmizp11NdN/Em2K4oEn1sFj/aXPyVxEckVkWPtj4HTgnb4ed6BZ45kT8lwnCxpCYkha/vZ6CCPkBd21aX7pWgoDdeQ74i8RRzLa52ekbxCUK+0eG65EfFYN0D83UvJ0BFwpyw2JKJmPA14WkbeAV4E/GWOeTsBxB4y/PPQ4u4JuLpHVGvk2aio9WXXNQDe9Py5eupbCwF7GZreE3CQErK5+vl4lNXunoEEz88rLwL7DUx499bPvarzZnbABSZkmep/3i4tvpTBQx+i8uj6fZ6DWyXb3GSOm//uW9PmMxpiPgaMTEMuA9a8tb3U8LgzsY+Z1OsfKUGLy/Ixr8DFh6jHdbjezNHrf996MwsyRZhrIwBVw4XN6O45wSelaVhfPZtKX/oeqd98F8vE7QudDjHa28LtDtQu+XWAkY3xeql1uRnsD7HFbZcAs114gej/xxPRXB0dYNp8Q2IvJMFT4uvYn7+tc4fHMwtndh0zJ0sURS+jJmmQLtGtij1YumQnOzvm4Z5beksJoVCrMXdjb7og9GzPlaJrfqybD1bUaZuqp3+SNF36P1yG4fSOZ8vVTO2OypynYtnUrv7/3ZuaU9r7Rd5yvBenhZkkXly8DoGrnTtZusDurhTVgJkuu8TI6UEeFYzgAM+wPzWjVGQcO97Nnj5tad//1ox8INJn3oDYokWs9uUq0H/10VtR1J5x0BiecFHlq3naTDz6YSyIkcpergvA7GUXjcDQz8+pbY9p23IQJMW0X03lz8ohlHjTBMLP05pj//352eRkbii+llhGM9EJNH9opg+9aFa+v/8/p/P31Z3p/8jjpRFvdWF98ScfjArSeXA0evb2lYDz8Dh+F/gbAqp7oeRRsqJlXL+f0k7+WjNASxi09TMLWjVPO6jq/jenjpGPd0WTejQpH51DfuR6tJ1cKOrtXGqcLk20lp1xXJWTGP3r0hG+eltDYEqH978sYYF0Pe6LJPApPiafzsVavqD4ydle9rIHaNSMOGXbvGofTzaxFN+HxeJhl38Q7mSbl1LNfXme9zIHD2yiUsFvaid2fP9OKMcvZ80jPDJ8bTGcqPPGrBzA5p5lFg2yuJU3m0dgNO7HMm61UTy5fuppxphq3e/BX1+Xn1VAYqGH2otAbeX/r7J+R43VR4KhIynl/teAmLryy86bbP7v8Wk749v9G3rZoJUdOHM5PryqNuD7YovIi9nd3Tv3wrenn8YsFy7vZY2DSBtAIgkvii8sG16ezGrhm90PptT/MXhj5PraHfnEqC5ZN7ddYvnTcSTzxh5cj3onofy+6vJs9Y/uKZPp4Y4wDxx3JsccdBZltPPLII306Vk80mYe57poZkGm12I/0JaeEoVR/KPTX0xQYnuowhrSfzba+OTQ1NSX9XJrMw7Rmdna9uqQ8ddOtqr677fpSnG4nc67o/HZ13aJrMI5MFpV7UhZXuCce/g3VVTu5aO6ChB53ZtlNUdd5PLeR43ABkScB8/t8mEAAp9vd64E4Oz79GIBJXzgw6jb/eWMTAQRHhJKy5DqgGdzuPV3W1dXWcsvK+7j0svMB8DqcNGVkcsxXv9xlWwPUZ+YwvLUzoYrDT0NtLVs/fJuKT3bj76FLemNGJn6Hk+EtncfId3hpxrrpuDfKHZf8PR04gSSZXWWimTZtmtm8eXO/n7cnNxbNpMFt9St3t+6g6Lo7k37O5YsuJXfESOYuKAHg2iIPAYfVK8BhWtg14kjG7PsvDl8je0ccSUNOHo9/6QsAjGz1U5MZ++24Dmis45Pc4RzYUMfHecPJ9fn5ccXjvJv7Zbbkj6XBaY1Ocwf8eB29v82XUiq6v375UL6Y17uRoCLyerQbAGkyDxJcV57oHiwP3HcXH295F3GO4pXDv8FrE/Trr1JD0Rn5bdxzbO8mlu0umWs1i221ZzbtI+YO3G9Mn49XvmAeOyedymN2KZpJ/2P9KKWGtN2tyak/12Ru220n8kyfm59deHGvjlF29VzuP/mX7Mt0wrcHz7zMSqn+8/i0Y5NyXE3mwLqSuSDWHc17c1Pmi295nEePPgDOjG3yf7dp5XJuIJ8a9jCGY9nMbgoYxR4EQwtZZNJKXU02n1ZPJse9hykHVFFVMYyPdkwlY6QPl2nAeNvw+zPw+11kZDTQuG8cX5zyOlt3Hoq/dhR+AfEbAv6R4A1wwplPsuWTQ6irHk1W7l5ysurYtW0Koyd/wu6PJ1N40Ae0VJ1Ak28Pfp8w8cC3qfjoGHD5QARn/m7yc3fT5suioW48i6++D4BNf3uB3//pbg760hbqW/Jx1h/HvoY9SFYLxpFJQIbhbNvH6Mmfs+edw5DsejLbxjJ28gQ++/wtTEBw5DbjzqonO6uOuvrxTD9+Hn//5x+ok8+gNYvFi+/i2hX/h2nLxDgCTJrwLrvrx9LUsB9nn3QBf3z5Fhx1I5g+/RccffxXKLvyEgJZrZSUr6N86UXgDIDXRdbYz5g88gf894PXyB82jpaszbTUTCArMIwWfwuLl1rbS0YbUpMP2X7MsH346/anKeswhrW9ghEfmaN3Mi6/ku1Vh3H0gdP54L1/c8jhx/Cfj/+Io3EU++1/KNv2/Asac3GSScaYD/nON1fyp8cepI7PGZ1dyBnTz2PKEUd0vC9uva6EhqYcisquBqB08Xxco3aCcRBoyaEgdzJtXi8jR4/l/AsuYdWKa5h/1XWsufNEXI4AVTsPg4xWcloncsWSG7ippJhWv4NF5R6Wlc3gjFPP48U3rqPJdwBfO/x0PtjyONtbDmCMqeJH5y3imT/ez2fVW5GAC2duI0ULfsPya+dx9aLbeOzBjWz56B8sLrqLG66bh5c6iq65h+f+8ihN9Q28s+0pLjxvZUdj6aZNz/Lux0+Qm1tNU2AiB+aexvsfVpI7bBfZmbnMmleOMW04HJl4PB5EAhx6UCGnf/vHrFlbxsETD8UhwtbqTQAcMPENJDCDc356MfffcxN19XuZPdea/KusxEPA4efwwly+9YsraWmoZ9KokdxQvorzfvY97ti4EbBun9fqyODAMSM4/PAjOPxLR/PXp/5Ibk4ulXWfsbtuEw5HK1de/Huef+lPbK/YQt2OLeQW7CErs5Gad7/KLqcLET+z51/M57s+4aPqNzn64JP47Yb7CQRcuDAceMjzfOmYq5gy+Rs4xIHf34aRbDbePYe21olknLws7hwTC60zp7N+fIzP2zE7XCxOe+gfvD0m+t3AF5qlHME7+HGSYfcaaK538+Z7X2Hxwnv6FLNSfbWuZB6VMpqRppJLlq7reYckaf//6207VdmScvxOH0dMyOLHMxZ2WV+6pJyA3Q99rK+NOeXXxnzs2xcvptpllXnH+n3scroixtr+3Ol3RbwDWVNTEzfccANZWVksXNg1xlhpnXk3NhRfAvYcLLEm8vk3P8JDxxwMERL5RnMumYROznP6tz4Mef6ds3sZrFIJNGvpbakOoV/8+Eff4cHHnog7kYeTAXubDMuQT+Y77URe4It+15R21y72sOqU78ExB4csv9HMpZDQAUavbDqRomu09K1UsjmN4AfcrsgTfR121FQ8R/XvyNRUGNLJfP2SK8Bpzcp/cVn3cziMf/5NOOV7ofubn5NHQ8fzgF847bStAJzyrQQHq1Sacvnd+JyRB93EYvIEP/WVdXz/Ak/igoqgu6FTJx8yiec/iH6f177eASkWQzqZV9iJfIzZhUjkOcfKF1/H6lO+HbJspZnD2KD5zRv2ZXD2D95LXqBKpbHxI1rw1u3t9f7/N6c4gdH0ziHHnsDzHzyU0hiGbDJf45kD9k2aL17a9U4tJmAofPEtCErkPzG/4bs83vH8k08L+NUvNyU9VqXS2a+uKEt1CGlhSCZzYwy77ESe56/ssv7eu9aw4MDQu4TcZ/63Y/6Iur1Z5LKMX/3ye132VUqpVBiSyXy9Zy7IGCTg5Mqy0C5ZpcU3sebkUzqef8P8jZlYU5dWV+Vxzrlv9WusSqmBYcrUCVS+Hd+t8frTkEvmNxRfhi8wAZeA07s9ZN2ld9zGg0GJ/A5zPjk08e/3pnDlxU/3d6hKqQHk8K+cxs4P7mTKsUd1XRlIfbfFIZfM842TCpeX8WYPs663SuXl1/8M53F+Hjz4qo7t7jc/5JXXvkLR1Q9ozxSlFADnX3NDn/ZP5iDNIZXM15VcTKVzDE6/i/2POImnHnwU99gFvH3cHF4Uq0Q+x9yC+8URnLL0I03iSg1xIqkvccdqyCTzaxfOoC3LuvHEONnDG1uf4rCvvclceRSAKf6POf0vf2PJil+DJnGlFPCNn5zFI49YVawud1aKo+nekEjmKxbOAtcXAC95Xic5h7+AmXAUc6Xz5hMXflrJL1asTl2QSqkB58gjj2dH3T62v/ICo8fuF3W7QArmuAqXkGQuImcCtwJO4E5jzPWJOG5f/On+h9nxwd+pkFGQNR5oo2DUVvYcmcE1rMcrGQD8bNM73LDwp3By+g/3VUrF78wTzoQTzoxx6+SP9Iymz8lcRJzA7cBpwA7gNRF50hizpa/HDueZdzlkN5CVNxx/RhZNOTl8XpBFW3YGzd7R5EkNzxacwGhTTW3hZHwTpgAwIrCXfY5RXY435+k7KV6upXGl1OCXiJL5ccBWY8zHACLyIHA2kPBk/u5ph/H3YT3fbmmPhN4pqD2R5/pbOe9fzzN5/HAuuGgOnKyJXCmVHhKRzCcCwR22dwBfCd9IRGYAMwD233//Xp3okM8rGTXqebJb28hqayWjpQ28PupGjiZrH+Saj2h1TsDdUI3D1Yxvr5evHPUdpv/fjzoPcmqX0JRSqm9MINURJCSZR6ok6tIaYIzZAGwA6+YUvTnRstmpn1BHKaXi1R+zJkaeKjA+O4DgZt5JwM4EHFcppVSMEpHMXwOmiMgBIpIBnAM8mYDjKqWUilGfq1mMMT4RmQv8Batr4kZjzLt9jkwppVTMEtLP3BjzZ+DPiTiWUkoNNsPHWFNq75+TE3F9MudkaTckRoAqpVQyDRs3nqvmzyN7WH7KYtBkrpRSCZA7anRKz5+IBlCllFLdGCxdE5VSSqWYJnOllEoDmsyVUioNaDJXSqk0oMlcKaX6STL7m2syV0qpNKDJXCml0oAmc6WUSgOazJVSKg1oMldKqTSgyVwppdKAJnOllEoDmsyVUirJdKItpZRSMdFkrpRSaUCTuVJKpQFN5koplQY0mSulVBrQZK6UUv3E7XYn7dh6Q2ellEqyzMxMTj31VA477LCknUOTuVJK9YMTTzwxqcfXahallEoDmsyVUioN9CmZi4hHRD4XkTftn+mJCkwppVTsElFnvtIYc2MCjqOUUqqXtJpFKaXSQCKS+VwR+Y+IbBSRkdE2EpEZIrJZRDZXV1cn4LRKKaXaiTGm+w1EngPGR1hVBLwC7AYMUAYUGmMu6Omk06ZNM5s3b44/WqWUGsJE5HVjzLRI63qsMzfGnBrjSe4A/hhnbEoppRKgTw2gIlJojKmwn34feCeW/V5//fXdIvJpL09bgPVtYKDRuOKjccVH44rPQI0L+hbbF6Kt6LGapTsich8wFauaZRswMyi5J4WIbI72NSOVNK74aFzx0bjiM1DjguTF1qeSuTHm/EQFopRSqve0a6JSSqWBwZjMN6Q6gCg0rvhoXPHRuOIzUOOCJMXWpzpzpZRSA8NgLJkrpZQKo8lcKaXSwKBK5iJypoj8V0S2isjCJJ9rPxF5XkTeE5F3ReQSe3nUmSJF5Bo7tv+KyBnJjFtEtonI23YMm+1lo0TkWRH50P490l4uIrLKPv9/ROTYoOP83N7+QxH5eR/iOTTomrwpInUicmmqrpc9vcQuEXknaFnCro+I/I99/bfa+0of4lohIu/b535cREbYyyeLSHPQtVvX0/mj/Y29jCthr52IHCAim+y4ficiGX2I63dBMW0TkTdTcL2i5YfUvceMMYPiB3ACHwEHAhnAW8ARSTxfIXCs/XgY8AFwBOABroyw/RF2TJnAAXaszmTFjdWvvyBs2Q3AQvvxQmC5/Xg68BQgwPHAJnv5KOBj+/dI+/HIBL1WlVgDHFJyvYBvAMcC7yTj+gCvAl+193kK+HYf4jodcNmPlwfFNTl4u7DjRDx/tL+xl3El7LUDHgLOsR+vA2b3Nq6w9TcBxSm4XtHyQ8reY4OpZH4csNUY87Expg14EDg7WSczxlQYY96wH9cD7wETu9nlbOBBY0yrMeYTYKsdc3/GfTZwj/34HuB7QcvvNZZXgBEiUgicATxrjNlrjKkBngXOTEAcpwAfGWO6G+Wb1OtljHkJ2BvhnH2+Pva64caYfxnrv+7eoGPFHZcx5hljjM9++gowqbtj9HD+aH9j3HF1I67Xzi5Rfgt4JJFx2cf9MfDb7o6RpOsVLT+k7D02mJL5RGB70PMddJ9cE0ZEJgPHAJvsRZFmiowWX7LiNsAzIvK6iMywl40z9ghc+/fYFMV2DqH/YAPhekHirs9E+3EyYrwAqxTW7gAR+beIvCgiXw+KN9r5o/2NvZWI1240sC/oAytR1+vrQJUx5sOgZf1+vcLyQ8reY4MpmUeqL0p6v0oRyQMeBS41xtQBa4GDsKYxqMD6mtddfMmK+2vGmGOBbwMXi8g3utm232Kz60K/CzxsLxoo16s78caSlBhFpAjwAffbiyqA/Y0xxwCXAw+IyPBknT+CRL12yYr3XEILDf1+vSLkh6ibRokhYddsMCXzHcB+Qc8nATuTeUIRcWO9UPcbYx4DMMZUGWP8xpgAcAfWV8vu4ktK3MaYnfbvXcDjdhxV9tez9q+Wu1IQ27eBN4wxVXZ8A+J62RJ1fXYQWhXS5xjthq/vAOfZX6uxqzH22I9fx6qPPqSH80f7G+OWwNduN1a1gitsea/Zx/oB8LugePv1ekXKD90cL/nvsVgq+wfCD9Y8Mh9jNbi0N658MYnnE6x6qlvClhcGPb4Mq+4Q4IuENgp9jNUglPC4gVxgWNDjf2LVda8gtPHlBvvx/yO08eVV09n48glWw8tI+/GoPsb2IPDLgXC9CGsQS+T1AV6zt21vnJreh7jOBLYAY8K2GwM47ccHAp/3dP5of2Mv40rYa4f1TS24AXROb+MKumYvpup6ET0/pOw9lpREmKwfrBbhD7A+cYuSfK4Tsb7W/Ad40/6ZDtwHvG0vfzLsDV9kx/ZfglqeEx23/UZ9y/55t/2YWHWTfwU+tH+3vykEuN0+/9vAtKBjXYDVgLWVoCTcy7hygD1AftCylFwvrK/fFYAXq5RzYSKvDzANa8rnj4DV2KOpexnXVqx60/b32Tp72x/ar+9bwBvAWT2dP9rf2Mu4Evba2e/ZV+2/9WEgs7dx2ct/DcwK27Y/r1e0/JCy95gO51dKqTQwmOrMlVJKRaHJXCml0oAmc6WUSgOazJVSKg1oMldKqTSgyVwppdKAJnOllEoD/x+SVUmdfRSCHgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "continue_train = True  # set to True to continue to train\n",
    "num_big_steps = 100     # number of small steps\n",
    "num_small_steps = 200  # number of big steps\n",
    "if continue_train:\n",
    "    for k in range(num_big_steps):\n",
    "        for j in range(num_small_steps):\n",
    "            dXY_list = np.append(dXY_list, np.zeros((rep, 1)), axis=1)\n",
    "            dX_list = np.append(dX_list, np.zeros((rep, 1)), axis=1)\n",
    "            dY_list = np.append(dY_list, np.zeros((rep, 1)), axis=1)\n",
    "            for i in range(rep):\n",
    "                minee_list[i].step()\n",
    "                dXY_list[i, -1], dX_list[i, -1], dY_list[i, -1] = minee_list[i].forward()\n",
    "        # To show intermediate works\n",
    "        for i in range(rep):\n",
    "            plt.plot(dXY_list[i, :],label='dXY')\n",
    "            plt.plot(dX_list[i, :],label='dX')\n",
    "            plt.plot(dY_list[i, :],label='dY')\n",
    "            plt.title('Plots of divergence estimates')\n",
    "        display.clear_output(wait=True)\n",
    "        display.display(plt.gcf())\n",
    "    display.clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save current results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current esults saved.\n"
     ]
    }
   ],
   "source": [
    "overwrite = False  # set to True to overwrite previously stored results\n",
    "if overwrite or not os.path.exists(chkpt_name):\n",
    "    minee_state_list = [minee_list[i].state_dict() for i in range(rep)]\n",
    "    torch.save({\n",
    "        'dXY_list': dXY_list,\n",
    "        'dX_list': dX_list,\n",
    "        'dY_list': dY_list,\n",
    "        'minee_state_list': minee_state_list\n",
    "    }, chkpt_name)\n",
    "    print('Current esults saved.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the ground truth mutual information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ground truth is 4.982193620464953 nats.\n"
     ]
    }
   ],
   "source": [
    "mi = g.ground_truth * d\n",
    "print('Ground truth is {} nats.'.format(mi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply moving average to smooth out the mutual information estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mi_ma_rate = 0.01            # rate of moving average\n",
    "mi_list = (dXY_list-dX_list-dY_list).copy()    # see also the estimate() member function of MINE\n",
    "for i in range(1,dXY_list.shape[1]):\n",
    "    mi_list[:,i] = (1-mi_ma_rate) * mi_list[:,i-1] + mi_ma_rate * mi_list[:,i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the mutual information estimate after different number of iterations. The red dashed line shows the ground truth, and the green dotted line is the number of iterations where 90% of the ground truth is reached."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hU1fbw8e9OCITQewuY0KWEAKGKEKV3RQSVqlfAwhULXhH1osJPUfG+ih0UUUBAAUWUjgLSIRA6oUuXXtOT/f5xJsNMMkmGZCZnZrI+zzMPZ05dMxlmzdn7nLWV1hohhBAiPT+zAxBCCOGZJEEIIYRwSBKEEEIIhyRBCCGEcEgShBBCCIcKmB2ArbJly+qQkBCzwxBCCK8RFRV1UWtdzh379qgEERISwrZt28wOQ/iIk9dOAlC1RFWTIxHCfZRSf7tr3x6VIIRwpUE/DwJg9dDV5gYihJeSBCF81uttXzc7BCG8miQI4bM6VO9gdghCeDWPTxBJSUmcOnWK+Ph4s0MRuRQYGEhwcDABAQF5cryjV44CUL1U9Tw5nhC+xuMTxKlTpyhWrBghISEopcwOR+SQ1ppLly5x6tQpQkND8+SYTyx8ApA+CCFyyuMTRHx8vCQHH6CUokyZMly4cCHPjvlW5Ft5diwhfJHHJwhAkoOPyOu/Y7uQdnl6PCF8jdxJLXxWzMUYYi7GmB2GEF5LEoQXWL16NT169MgwPzo6msWLF+don++88451+vjx4zRo0CDH8XmqEb+NYMRvI8wOQwivJQnCRZKTk/P8mFkliOzisU0Qvuqd9u/wTnvff53Cd2mtOXc11rTjS4Jwwvjx46lbty4dO3bk0UcfZdKkSQBERkYyduxY2rVrx8cff8zff/9N+/btCQsLo3379pw4cQKAoUOHMm/ePOv+ihYtChhnBpGRkfTt25e6desyYMAA0kb4W7p0KXXr1qVNmzYsWLAgQ0yJiYn897//Ze7cuYSHhzN37lzefPNNhg8fTqdOnRg8eDDTp09n5MiR1m169OjB6tWrGTNmDHFxcYSHhzNgwAAAUlJSGDZsGPXr16dTp07ExcW5583MQ62rtqZ11dZmhyFEjq3cdZohn/zJ7hOXTTm+V3RS24mMzDivXz945hmIjYVu3TIuHzrUeFy8CH372i9bvTrLw23bto358+ezY8cOkpOTadKkCU2bNrUuv3r1KmvWrAGgZ8+eDB48mCFDhjBt2jSee+45fvnllyz3v2PHDvbu3UvlypW55557WL9+PREREQwbNow//viDmjVr0r9//wzbFSxYkLfffptt27bx6aefAvDmm28SFRXFunXrKFy4MNOnT3d4zIkTJ/Lpp58SHR0NGE1Mhw4dYvbs2UydOpV+/foxf/58Bg4cmGXsnm7P+T0ANCjve81nIn/YfeISAAdOX6FhtdJ2y4Z9sYYTF2+69fhyBpGNdevW0bt3bwoXLkyxYsXo2bOn3XLbL++NGzfy2GOPATBo0CDWrVuX7f6bN29OcHAwfn5+hIeHc/z4cQ4cOEBoaCi1atVCKXVHX9S9evWicOHCTq+fJjQ0lPDwcACaNm3K8ePH73gfnmbk4pGMXDwy+xWF8FCbDp4H4OuVB+zm/3M11u3JAbzxDCKrX/xBQVkvL1s22zOG9NKafDJTpEiRTJelXdZZoEABUlNTrftLTEy0rlOoUCHrtL+/v7XvIKeXhNrGY3tcIMu70dPH4QtNTB90/MDsEITIMa0112KN74oKJex/9P0WdSJPYpAziGy0adOGRYsWER8fz82bN/n9998zXbd169bMmTMHgFmzZtGmTRvAKGMeFRUFwMKFC0lKSsrymHXr1uXYsWMcOXIEgNmzZztcr1ixYty4cSPT/YSEhBAdHU1qaionT55ky5Yt1mUBAQHZxuHtmlVpRrMqzcwOQ4gciU24faFJ23qV7Jbl1S1FkiCy0axZM3r16kWjRo3o06cPERERlChRwuG6kydP5ttvvyUsLIwZM2bw8ccfAzBs2DDWrFlD8+bN2bx5c5ZnHWDULJoyZQrdu3enTZs23HXXXQ7Xu++++9i3b5+1kzq9e+65h9DQUBo2bMjo0aNp0qSJddnw4cMJCwuzdlL7ouhz0USfizY7DCFy5Eb87R9wAf72X9Vz1x/JkxhUdk0oeSkiIkKnHzBo//793H333SZFZLh58yZFixYlNjaWtm3bMmXKFLsvW+G8vPx7Rk6PBKQWk/BOR85d45mpRj9m8cIB/DS6k3VZ1wm/k2r56l7+3x5RWusId8TgfX0QJhg+fDj79u0jPj6eIUOGSHLwEh91+cjsEITIsauxt/sqr8fZNwd3CAtm+c5Tbo9BEoQTfvjhB7NDEDkQXjHc7BCEyLFLN+wvKtFao5Ti7ws38iQ5gPRBCB+29fRWtp7eanYYQmQr5sxVLl63TwjpnycmG1ckTvw57/rV5AxC+KyXV7wMSB+E8GybDv7DuLlG3+v/PdacsLtKU7CAPxfTnUEkJKdQKMCfo/9cz7PYJEEIn/Vpt0/NDkHkc1dvJVA8qCB+WVyXmpYcAF77wbgU/dE2NTOcQSRZziD8/RQpqXlzcZEkCOGzpMSGMMuF63FcuhHPqGkbKFY4gHk2VyA5Y/a6w9SoUNxuXloTU14lB5A+CKc88cQTlC9fPkNJ7MuXL9OxY0dq1apFx44duXLlCmB0Jj333HPUrFmTsLAwtm/fDkBMTAxNmzalUaNGbNy4ETCqrnbo0IHYWPMqNmYnfdE/Z4SEhHDx4kU3ReScDSc3sOHkBlNjEPnTwI//YNQ047N3I90VSMkpqfy8+RgpqansOJb5/5EjlqakooHGGO4JSSluijZzkiCcMHToUJYuXZph/sSJE2nfvj2HDh2iffv2TJw4EYAlS5Zw6NAhDh06xJQpU3j66acB+Oqrr5g4cSLz5s2zVoT94osvGDRoEEFBQXn3gvKJsavGMnbVWLPDEMLOgs3H+HL5Pn6LOsGYmZuzXT/tjuoRX61l9d4z1vkfP+H+SsWSIJzQtm1bSpcunWH+woULGTJkCABDhgyxVm5duHAhgwcPRilFy5YtuXr1KmfPniUgIIC4uDhiY2MJCAjg6tWrLFq0iMGDB2d67CNHjtClSxeaNm3Kvffey4EDRtGu3r178/333wNG4km7I/rw4cN06NCBRo0a0aRJE2u5jg8++IBmzZoRFhbGuHHjrPufOXMmzZs3Jzw8nBEjRpCSYvxK+fbbb6lduzbt2rVj/fr11vUvXLjAQw89RLNmzWjWrJl12aVLl+jUqRONGzdmxIgR2dawygtf9fiKr3p8ZXYYIp+ITUjm0Nlr2a6XdkYRl+DcGDKpNv+X3l2wgyKFClC8cAB1q5SiViXHVR1cxesSROT0SKZHTwcgKSWJyOmRzNw1E4DYpFgip0cyd49RduJa/DUip0eyYL8xnsLF2ItETo9kUcwiAM7dPJerWP755x8qVTJqpFSqVInz543Ki6dPn6Zq1arW9YKDgzl9+jTPPvss//vf/3jqqacYO3Ysb7/9Nq+99lqWhfmGDx/OJ598QlRUFJMmTeKZZ54BYMqUKbz99tv89ddffPjhh3zyyScADBgwgGeffZadO3eyYcMGKlWqxPLlyzl06BBbtmwhOjqaqKgo1q5dy/79+5k7dy7r168nOjoaf39/Zs2axdmzZxk3bhzr169nxYoV7Nu3zxrPqFGjeOGFF9i6dSvz58/nySefBOCtt96iTZs27Nixg169elnHwjBTnbJ1qFO2jtlhiHziwfeXMfLrdRk6l9NLsRTQ/PZP54bDHRJZ2+55iSIFaVK9HACfPtkmB5E6Tzqp3cDRr2elFNWqVWO1pZrs4cOHOXPmDHXr1mXQoEEkJiYyfvx4ate+/WG4efMmGzZs4OGHH7bOS0hIAKBChQq8/fbb3Hffffz888+ULl2aGzducPr0aR588EHAqOkEsHz5cpYvX07jxo2t+z106BC7du0iKiqKZs2MgnZxcXGUL1+ezZs3ExkZSblyxoewf//+HDx4EICVK1faJYzr169z48YN1q5dax3YqHv37pQqVSr3b2QurTlujNPRLqSdyZGI/OTlGRszzOs8/nd+eqkjxYMKMn/TsTvaX53KJe2e34pPJqhQ3nx1e12CsL2mPcA/wO55UECQ3fMSgSXsnpcNKmv3vGLRirmKpUKFCpw9e5ZKlSpx9uxZypcvDxhnDCdPnrSud+rUKSpXrmy37WuvvcaECROYPHkyAwYMICQkhLfeeotZs2ZZ10lNTaVkyZLWgX3S2717N2XKlOHMGaNdMrNmHa01r776KiNG2I/P/MknnzBkyBDeffddu/m//PJLpmc1qampbNy40eGYEzktUe4u41YbTWlyH4TIS2cuO77g5KXvNjL16Tv/sVIk0P5r+lpsImcu38pRbHfK65qYPEmvXr347rvvAPjuu+/o3bu3df7333+P1ppNmzZRokQJa1MUwJo1a6hSpQq1atUiNjYWPz8//P39M1zJVLx4cUJDQ/npp58A44t+586dAGzZsoUlS5awY8cOJk2axLFjxyhevDjBwcHWvpCEhARiY2Pp3Lkz06ZN4+ZNY4CR06dPc/78edq3b8+8efOsTWOXL1/m77//pkWLFqxevZpLly6RlJRkPT5Ap06drCPYAdbk1bZtW2tyW7JkifWKLjNN6z2Nab2nmR2G8FC3EpJcNt7z5ZtZNysBnLh4k+jjmV+1VK1sUYfzQ8oXzzAv+vgl54PLBUkQTnj00Udp1aoVMTExBAcH88033wAwZswYVqxYQa1atVixYgVjxowBoFu3blSvXp2aNWsybNgwPv/8c+u+tNZMmDCBN954AzD6GMaMGcNDDz3E6NGjMxx71qxZfPPNNzRq1Ij69euzcOFCEhISGDZsGNOmTaNy5cp8+OGHPPHEE2itmTFjBpMnTyYsLIzWrVtz7tw5OnXqxGOPPUarVq1o2LAhffv25caNG9SrV48JEybQqVMnwsLC6Nixo/WM6M0336RVq1Z06NDBrjjh5MmT2bZtG2FhYdSrV48vv/wSgHHjxrF27VqaNGnC8uXLqVatmtv+Hs6qXqo61UtVNzsM4YGSU1Lp8/5yhnzyJ9uPOnc59u9Rf/PKjE0Ol42dtcXh/PRemeH4qqXI+pX5YHBLhnXIWOk4MMDfqX27g5T7FnkqL/+eK4+uBKBD9Q55cjzhPWavO8x0m07iZW90z3abzuONwcJmjrqfcsULO1yWUy/2DKNzeFWu3kqg//9W2i1b9kb3DPtPWx9AKeW2ct9yBiF81oS1E5iwdoLZYQgPdPH67SF1SxctlMWaGQ38+I9s1/nppY5OJZ00e09eBqBkkUKM6t6QqmXsBxVr37CK3fPMmqNczes6qYVw1owHZ5gdgvBQB23uV3iwRajL9188qOAdrW8zdDzdmlTjrnJFeXH6RtIu+/jPA+Gs2n3auk76EebcxSvOIDypGUzkXF7/HauWqErVElWzX1HkO6Hli1mn0z6XKampTPp1p1M3u+0+cTnTZYte7XLH8XRvat9nl1ZuqV5Vx5eLF5AEYQgMDOTSpUuSJLyc1ppLly5Z783IC0sPL2Xp4YwlUoRYFn17wJ20KqlnLseyYucpvli2N9vtR393+16HCzbNVQ+3qk7BAnfWqdytSTXqVCmZ7XqvPHB7AKwC/nlzSbnHNzEFBwdz6tQpLly4YHYoIpcCAwMJDg7Os+NNXGfUxupS885/0QnfdSvhdvE8fz9FYoqRIG7GG/Ov2Qz1mWbTwX8c7uvyzXi7Pol/ta97x/GM6FQvy3LgaZpUL2udLuCXN7/tPT5BBAQEEBrq+jZC4fvm9J1jdgjCA83+67B1OsDfjyRLgkirrHrqkv1NaPM3HWXKiv0Z9nPiwo0MyST9zaIdw4JZsSvr4UEDHJwNVC9fjAB/PwbcW+v2egVuJ4W8amJya4JQSr0APAloYDfwuNY6+ztKhHCB3N4pL3xfwQJ+1iam8iXsL12NT0ym93vLMt122Jdrs70z+qVeYdYEUbigP3GJt0t2v/pgYxJTUvB3cDZQJDCA38Z2TRfr7aarvGpiclsaUkpVAZ4DIrTWDQB/4BF3HU+I9BbFLLIWZhQCjOalnzYetT4vXKiAtZz2Bwt3Wud/+8cBjl+4me3+zl25fSf2J/+6J8Ny2zOKn0Z3IqSc0TnesFppIhtUplMj5y+iKOB3e195dRWTu5uYCgCFlVJJQBBwJpv1hXCZDzd+CEDPOj1NjkR4ipU77Zt7ggreThC25qw/wpz1R7Ld3xtztlqnAws6/jqd+nQ7ChXwI8Dfj8+GteGHvw7Tr/Wd3+Fvm2wK51GxPrelIa31aWAScAI4C1zTWi9Pv55SarhSaptSapt0RAtXmtdvHvP6zTM7DOFBPl92uxLxR4+3JjDAn40H/6Hfhyuc2r5l7QqZLrsVn+RwfrWyRalQ0hgQrIC/H4Mja2eaTLKTdlOfM53aruDOJqZSQG8gFKgMFFFKDUy/ntZ6itY6QmsdkVZeWghXKBtUlrJBZbNfUfica7GJxCVmPSBP0cAA9p++al0/K2l3LlcpnfnIjyE291a4y+R/3cP4R5q5/Thp3Hme0gE4prW+AKCUWgC0Bma68ZhCWKUNFNXn7j4mRyLyWtoZQVblLsoVd/6enLC7SjP16XbEJ6Xw69a/rVc+2Sqcw7OCO1GueOEMdaDcyZ09HSeAlkqpIGU0nrUHMl4rJoSbTN48mcmbJ5sdhjDR9D9jmP5nDKlac/VWgnX+k+3rZtrM80KPhhnm+Vk6iAMD/PltbNcMTTzvDWzhwqg9h9tSntZ6s1JqHrAdSAZ2AFPcdTwh0lv4yEKzQxAudvVWAp8v3cuoHg0pUijA4Tq2VRdmrzPuebgWm8iDzUOs8/u0dNxJ/GCLULo0rsb/+2233fz0CSE1XWWH8FDfbMp067VSWutxWuu6WusGWutBWuuE7LcSwjVKBJagRKB7B3UXrrNq1ynr3cyZ+XrlAdbsO8uL32Yc1jPNP9fiMsxbvP0EKam3v9T9LWcE/VvXsFuvexOjJtJ3I++zm5/VaIkP2CQeX+PxtZiEyKm5e+Yyd89cs8MQTth+9CLvL9zJ/xbtynK9/aeMkQqPX7iR5b4cuXTT+H1qW0r7ifZ16RJ++16EQpbBeSqWsu+MLlUk85LgzWqWzzJmb+bxpTaEyKkvtn0BQP8G/U2ORGTlVkISr84yRlpbf+BcpuulpGpO2YzFPPaHLdSpXIIhkXXs1jt45qrD7T9ZbDQbPd2lvt182/EgCtqUs/jppY5cuB7H3xdu0q5+Jbttvhh+L09P+QuAuk4U2vNWkiCEz1o8YLHZIQgn9Hk/w+1RDm0/an+fVNSRC0QduZAhQVQpbT/YTppzV42mp0s3Mq/2Y3slUvGgghQPKkiNihmbKatXuD1OtG1S8TWSIITPCgrI/Jp14R3mrDtMRI1y1KxUgtdnb3W4ju1QoO8u2MHek1ey3GeJdIP5hNp82RfKwfjPeVX2wgySIITPmrnLuOVmYFiG+zOFB0s//vK3f8aw4OVO2W7nzFCgAC1q2d8Nfe/duSvqmFUHtreTBCF81tfbvwYkQfiCPh841wyVE778BZ9bkiCEz1oxyLn6OsI8zgzvmVNPd67HFza1l7LyVv+ITGspZSaoYAFisynn4e0kQQifFeDv+EYqkXfOXoll9rpD/LtbQ4dt9ZsPnb/jff78n848+H7m4zSkaVmrAgH+ftzXoApXbiXwxGerM183iyJ8mfnhhfYOK8H6EkkQwmdNj54OwNDwoabGkZ8kJqcQl5hi7Qh+4dsNXLmVQJligRmuNgKYseagdbpSqSDO2oyvkJkgJ0tdKwXdm951R9vcicIFC+RJ/SUz+W73u8j3pkdPtyYJkTfenLvNrnT2FUv9ox9shvlMY1sS44HmIUwckH09o393a+BUHMM63J1hhLj/9G7EqO4Z6yyJzPl2+hP52uqhq80OId+JstzFnKo1Jy86HpHtly3HWLX7NI/cU9M674HmoVQomX2V0h6WM4LPh93LjfhEXpmxOcM69auWom+rjLWW2ocFO/UaxG1yBiGEcAnbMRXOXo5l+JdrM6wTc+YqXyzbx8Ez1+yakyqVCrK7muit/hF89HjrTI9Vo2JxwkNuF8gb1uFu6/Qznes72kTkgCQI4bOmRk1latRUs8PIN37465B1+onPVztc57lv1lunHRXme6NvE57tUp+WtStwd3Ap/tW+rnXZiI53Z1g/TZfGt+sp1awkBRpdRZqYhM+au9co1Des6TCTI8kfbKulOpKQlGL3PK0U98dP3GOd1+Zu+5pH/VrXoEzRQry/cCddLZVWbc0cdT8Hz1yjaGAAEx5tlqHfQeSOJAjhs1YOXml2CPlKdiO0/XeO41IZ2W3XPiw40/4D2xHWfLmqqlkkQQghcu2bVQf4ccORLNeJPn7J4fyc1D8SeUP6IITP+nzr53y+9XOzw/BJcYnJ1nGZtdZZJofB7WpnuS9JEJ5LEoTwWYsOLmLRwUVmh+FzTl68yQPvLaPHO0tI1dru6iVbHz9xD8ve6M4fe05nuq+eEXf5dDVUbyd/GeGzlgxYwpIBS8wOw6tprZn2xwG7S1Kf/GKNdbrrhMX0/5/jvp5alquJ0l92+tpDTazTtSvLFUeeTPoghBCZOnMllrnrj7Ax5h+e696QBZuOZrl+j6bVaFm7Aicu3rSO+1yjYnG7ddrWq0T1Cu34fNm+DKW3hWeRBCF81sebPgZgVMtRJkfieX7acIQWtcpTrVyxLNfbc+IyACcu3mT0dxuzXPfu4JI827UBfkrZXVHkqA5ScJmivPNY8xxELvKSNDEJn7Xq2CpWHVtldhgeJy4xma9XHWCYgzud03Om/EWaamWL4udgbIWCBfx52FL6Iv1obsKzyRmE8Fm/Pvqr2SF4pNkOCudlJj4xJfuVLOpWKZXpsgKWjuiHW2eskSQ8l5xBCJGPbDr4D3OzuV8BjJpJv0X9ba3Gmt64h5tmmNembuZDdz7cujo9I+6ip6XYnvAOcgYhfNakDZMAGN16tMmReI5xc7fZPZ+19hAD2tbKsJ5tzSRHggIzfnUUz6L5qEihAEZ2da5Ut/AckiCEz9p4KutOVV93/locUUcv0LWxUcPI0fCe3685yPc2g/a8O6AFDaplbCoa168pf+w+zfM9wjh16RY1KhYntHwxht5XJ0PSEb5DEoTwWfP7zTc7BFO99sMWTly8SZ3KJaleoTj//npdttu8Oivj+AoAretUpHUdowmpbpWSAHw5om2GAnzCt0gfhBBebNDkP/h86d4M86/HJnLCMmDP01P+AqC8zRVJT3Wq55LjFwrwp2n1sjzUMtQl+xOeRRKE8FkT101k4rqJZofhVuevxbFw63G7eZduxPOwzbCfYIzwdv5qHACLX+vKgy1c94X+zoAWDO/omoQjPIs0MQmfFX0u2uwQ3Co+XfPO8fM3+Gv/WQ476GsYNW09aaM1+PsZvwuXvdGdzuN/z/Y47w9qmetYhXeSBCF81py+c8wOwa22H71gnd578jIvTs+8U/7gmYxJA+C7kfdx9J/rvPVTVIZlFUsW5tzVOIo4uBNa5A/ylxfCS0Ufuz2+QlbJISsVSwVRsVQQrzwQTkqqZtKvOzOsU7igfE3kV9IHIXzW+DXjGb9mvNlh5Ng3qw6wM5NBdgCqlA7KcvvFr3Wlbb1KWa6T5v6GVejYKJhvnmkHwHsDW9DdclNbsaAAJyMWvibbBKGUqq2UWqWU2mN5HqaUet39oQmROzGXYoi5FGN2GDmSkmoMwvOfGZvs5icmp/D3hRsArNp9JtPt76lTAX8/P8b2aWw3/6sRbbM8bnCZoix7ozvhoWXp26o6P/+nM8ULS/2k/MqZc8epwMvAVwBa611KqR+ACe4MTIjcmtlnptkh5Nit+CSH83u+uxSA9g2rEHPmaqbb/7dfBABKKSLrV2b1XiOZhJTPunqrLT+lHFZiFfmHM01MQVrrLenmJTuzc6VUSaXUPKXUAaXUfqVUqzsPUYj8x/au5582HuF6XCLPT7td/mLVbsejtPVpEcq7A1rYzXvs3pqULlqI2S+0d0+wwmc58/PgolKqBhhXySml+gJnndz/x8BSrXVfpVRBIOtGUyFc6L9//heAt+972+RI7tzYH27/Jvt65QG+Xnkg03VfeSCc936Jpm+r6gzrcHeG5XeVK8bsFzq4JU7h25xJEM8CU4C6SqnTwDFgQHYbKaWKA22BoQBa60TA8eC1QrjByesnzQ7BaXtPXiYwoAA1Khbnj0zODhwpVzyQ+xtWoXXdigQG+LsxQpEfOZMgtNa6g1KqCOCntb6hlHLmNszqwAXgW6VUIyAKGKW1vmW7klJqODAcoFq1ancWvRBZ+Lb3t2aH4LS0y1SXvdGd935x7ga/nhF3MeBeoxKrJAfhDs70QcwH0Frf0lrfsMyb58R2BYAmwBda68bALWBM+pW01lO01hFa64hy5co5GbYQvkNrnf1KFoMspbkbh5ZlZNcGlCpayF1hCZH5GYRSqi5QHyihlOpjs6g4EOjEvk8Bp7TWaeUh5+EgQQjhLq+ufBWAdzu8a3IkWdt+9KJ1On3pi6c61ePL5fsAeLFnGJ3Dq9IopAw1K5XI0xhF/pRVE1MdoAdQEuhpM/8GMCy7HWutzymlTiql6mitY4D2wL7cBCvEnbgUl/lNZmaJTUjmwfeXAfDOgOaEh5S165C2Vbl0EA80D+H05Vu0rF2BiBrGGXbDu8rkWbwif1PZnd4qpVpprXN0H79SKhz4GigIHAUe11pfyWz9iIgIvW2bDD4ifJczxfHSLHm9G35KuTEa4QuUUlFa6wh37NuZTuodSqlnMZqbrE1LWusnsttQax0NuCVwIbxJSqqm2/8tvqNtJDkIsznTST0DqAh0BtYAwRjNTEJ4tNHLRzN6ufnjUV+9lXDHyWHqU1mXxBAiLziTIGpqrd8AbmmtvwO6Aw3dG5YQuReXFEdcUpzZYfD4p6vvaP3PnmxDtXLOl8QQwl2cSRBpRWGuKqUaACWAELdFJISLfNb9Mz7r/pnZYXBX+aJ2z4Nsymd/9HjrDBdK8TcAAB9NSURBVOvLFUrCUzjTBzFFKVUKeAP4FSgK/NetUQnhpeISk1m77yztG1ahgL8fWmv2n7IvqlejYnHGPNiYhOQUqpQuwrh+TSlfvDAz1x5iy+HzJkUuREbZJgit9deWyTUYd0cL4RWeX/o8AB91+SjPjvnr1uNM+yOGwAB/GoeWtRsbeuao+/n31+v5YHBLlE0HdOs6FQF4s79czyE8S7YJQilVEhiM0axkXV9r/Zz7whLCs2y1/LJvVrN8luvFWIb23H7sIu8s2GG3rFzxwsx5UYrmCe/hTBPTYmATsBtIdW84QrhOTs8cUlJTmf7nQR5sEULposaV3a/P3goYtZIyc/ZKLOsPnAOMCqpCeDtnEkSg1vpFt0cihIfYduQCP244wvlrcbzapzFbDt3uF/gt6m+a1yzP5ZvxBPj7U6NiceuyoZ/+aZ3+arkUDRDez5kEMUMpNQz4DUhIm6m1vuy2qIRwgWd/fxbgjq5kStWaNXuN4U5iE5OZ/mcMs9cdti7/ZPEeu/XTziiux0kle+F7nEkQicAHwGtYBg2y/Csd1sKjFQ4ofMfbvPbDFmvxvC2HztudPThy/locgyb/ke1+H2zhTIV8ITyLM7WYjgAttNYXs1zRBaQWkzDLpRvxfLPqQKZDeTrr/UEt+c+MTdbnFUoU5rt/32d31ZIQruTOWkzO3Ci3F4h1x8GF8BRPfLY6y+Sw9PVuTu3nyD/X7Z4PaFtLkoPwWs4kiBQgWin1lVJqctrD3YEJkVvDFw1n+KLh2a53/PwN4pNS7Ob9q31d6/Tkf92DUorRvRpRt0pJHrmnRqb76tG0GiM6GuNCf/tsJJ3Dq+YweiHM50wfxC+WhxBepUzhMpy9Gsvek5epX7V0huUrdp7izJVbXL1l38HcpHpZ+raqzjerDgBQp3JJADo2CqZjo2AAht5Xhy4TbhfgKxoYQGT9ShQs4E+fltXp01K66IT3y7YPIi9JH4RwpaSUVHq8swS4fbXRuv1n+WzpXi7fTMh0u0+fbEMtJ+oh/bnnNKcvxzLg3prSjCRMY8p4EEqpH7XW/ZRSu7l99ZKV1jrM5dHExEBkpP28fv3gmWcgNha6OWgHHjrUeFy8CH37Zlz+9NPQvz+cPAmDBmVc/tJL0LOncewRIzIuf/116NABoqPh+eczLn/nHWjdGjZsgLFjMy7/6CMID4eVK2HChIzLv/oK6tSBRYvgww8zLp8xA6pWhblz4YsvMi6fNw/KloXp041HeosXQ1AQfP45/PhjxuWrVxv/TpoEv/1mv6xwYVhifMEyfjysWmW/vEwZmD/fmH71VdiYblyp4GCYOdOYfv554z20Vbs2TJliTA8fDgcP2i8PDzfeP4CBA+HUKfvlrVrBu5bhRB96CC7ZjyDXo/3L1mndtSsqLo7xNvPSG3xkHfH+AdR69ANjRjafvfvksyefPXD42aN9e3jjDWO6a1eIS1dVuEcPGG0pRZ/+Ow/u7HvPjbJqYhpl+beHWyMQwg00sCfV+A/ewO95urQeybJVH2S5zYDjORo4UQif5cxlru9prV/Jbp4rSBOTcJU35mzlh5j/AVDTb2C26099qq2MwSC8ktlDjnYE0ieDrg7mCeExLl2Pp6bfQN4b2IJXZm7OsHzS4JYUCQwgtHwx6T8QIhNZ9UE8DTwD1FBK7bJZVAxY7+7AhMiNtPsRwkPL0qxmObYevmBdtuDlThQJDDArNCG8Rlb3QfwA9AQWWv5NezTVWmd/zi6ESUZ+vQ6A3amTGLhgIK1qV7BbLslBCOdkmiC01te01seB14FzWuu/gVBgoGWMCCFMlZicQmJyCnGJydZ5SSmpHDprjMkQpKpQp0wd2ocZ9y6ULR7Id/++z5RYhfBGznRSRwMRGAMGLcMYdrSO1tq52gN3QDqpxZ3oPP536/So7g3p1qSa3bwhkbV57N5aAKR9zqW/Qfgas2sxpWqtk4E+wEda6xeASu4IRghH9p26Qufxv3MjLonj529wMz6JI+fsax59/Ptuu+TQsFppa3IAIzFIchDizjhzFVOSUupRjGFHe1rmSSOucKmEpBTGz4tiSGQdShUpRNnigdZlL3y7AYC+k5Y7vb/3BrXgkXmPADCn7xzXBitEPuFMgngceAr4P631MaVUKDDTvWGJ/KbXxKUA1quN0kpjHLb0J2SmSukinL58y27eo21q4u/nR3jFcDdEKkT+kW0Tk9Z6H8Y9D9stz49prSe6OzCRv3Ue/zsnLtzgWcsVSZn5YHBLOoRVsZs39L46AIxpM4Yxbca4LUYhfF22ZxBKqZ7AJKAgEKqUCgfe1lr3cndwIn9IX2o7zbAv12a6zbNd6lOtXFHKFAvk5d7hvNzbOFvwpOKTQng7Z5qY3gSaA6sBtNbRlmYmIXLs4vV4Chbwo1jhAB6wNC9l5d0BLXh11mYm/+sea/ltR2w7oh/68SEA5vebn/uAhciHnEkQyVrra+muAJGfaSLHtNYM+HhVhvmTBrekSpkilC4ayKpdp3h/4U7rsibVy1r7JZzVKrhVrmMVIj9zJkHsUUo9BvgrpWoBzwEb3BuW8FXjf4pi3YFzDpc1vKuMdbp9WDDhoWX5euV+RvdulKNjjW49OkfbCSEMztwH8W+gPpCAUX7jGuCgOL0Q2cssOXzzTLsM88oUC+SVBxvj7+fMx1QI4WrZnkForWOB1ywPIe7YlkPnuRabaB2uM803z7QjuExRtx2312zjOopfH/3VbccQwpc508QkRI6duxLLG3O2AjDpV6NPoXezEIZE1nZ70bz2oe3dun8hfJ3bE4RSyh/YBpzWWsvodPmE1pofNxxl2h8HMiy7fDM+Tyqqjmo5KvuVhBCZyosziFHAfqB4HhxLmGzh1uN8vnRvluu83rdpHkUjhMiNrAYM+oQsLmfVWj+X3c6VUsFAd+D/gBdzEqDwHrEJyZkmh7f6R7Dj2EWe7lw/z+LpOqsrAEsGLMmzYwrhS7I6g3BF3e2PgP9gjEInfNjrs7fYjdqWZvK/7iEpOZUG1UrTMt3APe7Ws3bP7FcSQmQq0wShtf4uNztWSvUAzmuto5RSkVmsNxwYDlCtWrXcHFLkkbSy2iM61eOnDUe4fDPBbvlb/SPyPBk48kyzZ8wOQQivlumAQUqpLK8NzK4Wk1LqXWAQkAwEYvRBLMhquFIZMMjzPf7Zn5y5HJvp8odbVefJDnfnYURC5G/uHDAoqyamVsBJYDawGbij0Va01q8CrwJYziBGy1jW3i+z5FAiqCAzR91PwQL+eRxR5jp83wGAlYNXmhyJEN4pqwRREegIPAo8BvwOzNZaZ32JivBZsQnJDucverWLRyWGNP3r9zc7BCG8WrZjUgMopQphJIoPMEp9f+KOYKSJyfPEJiSz5fB5Cvgpxs/bDsDoXo24ciuBb1Yd8NjkIER+YVYTU1pi6I6RHEKAycACdwQiPE9icgoPvr8sw/x29StRsIA//VrXMCEqIUReyeo+iO+ABsAS4C2t9Z48i0p4hJ7vZhyn4fexXSng7x3F8yKnRwKweuhqU+MQwltldQYxCLgF1AaesxkPQgFaay13Rvuo4V+u4e8LNzPMn/9yJ69JDgBDw4eaHYIQXi2r+yC855tA5Nqt+CQ+WbKHf3drYJccXuwZRufwqiZGlnOSIITIHanmKjhw+gqjphljQP2554zdMm9NDgBJKUkABPi7vzCgEL5IEkQ+9/PmY3y5fF+G+d8+G0nl0kVMiMh1Os7oCEgfhBA5JQkiH0vV2mFyALw+OQA82eRJs0MQwqtJgsinrscl8vCkFdbny97oTnxiMit2neb+hpVNjMx1BobJjftC5IYkiHwoMTnFLjn8+JLRFBNYsAA9I+4yKyyXi00yyoIEBQSZHIkQ3kkSRD4zbu42Nh38x/q8YAE/SgQVNDEi9+k2qxsgfRBC5JQkiHxk08F/7JLDvNGdKFbYd6/weTriabNDEMKrSYLIBxKSUth25AJv/xRlnffl8Ht9OjkA9G8gxfqEyA1JED4uJVXTa6J9yYyZo+6nXPHCJkWUd67FXwOgRGAJkyMRwjtJgvBhZy7f4vHPVtvN+2L4vfkiOQD0ntMbkD4IIXJKEoQPSknVbIg5xwRLee40b/WPoHqF/FNC67kWz5kdghBeTRKED7kem8jDH67IMH/qU22pVq6YCRGZq8/dfcwOQQivJgnChzhKDvl5QJ+LsRcBKBtU1uRIhPBOkiB8RJfxv1unSxctxP891jxfNSc50vfHvoD0QQiRU5IgfMChs9dIGzh2UNtaDGxX29R4PMVLrV4yOwQhvJokCC+XmJzCyK/XAfBmvwha1algckSeo2ednmaHIIRXkwThpeITk+n9nv140ZIc7J27eQ6AikUrmhyJEN5JEoQXOn7+BiO+Wms3b/K/7jEpGs/1yLxHAOmDECKnJEF4mc42ndEAz3SpT6+Iu7AZM1xYjGkzxuwQhPBqkiC8yBtztto9/3L4vYTm8yuVstKlZhezQxDCq0mC8BLf/RnDlkPnAfjh+faUKRZockSe7+S1kwBULeG942oLYSZJEB5Oa82TX6zh1KVbAITdVVqSg5MG/TwIkD4IIXJKEoQH01rTZcJiu3kfDG5lUjTe5/W2r5sdghBeTRKEB5v4c7R1OrJ+ZZ7tWt/EaLxPh+odzA5BCK8mCcJDXb2VwOq9Z4D8XU8pN45eOQpA9VLVTY5ECO8kCcIDfbpkD4u2/Q3AA81DJDnk0BMLnwCkD0KInJIE4WHS3+fwdGdpVsqptyLfMjsEIbyaJAgP8sni3dbp/DIsqDu1C2lndghCeDVJEB5i25EL/BZ1AoAfX+pIiaCCJkfk/WIuxgBQp2wdkyMRwjtJgvAAZy7f4rUftgAw4dFmkhxcZMRvIwDpgxAipyRBmCw+KYXHP1sNwMC2tWhWs7y5AfmQd9q/Y3YIQng1tyUIpVRV4HugIpAKTNFaf+yu43mjlNRUnp3yFwD1gksxSAb6canWVVubHYIQXs2dZxDJwEta6+1KqWJAlFJqhdZ6nxuP6TXe+3kHf+wx7nOIrF+ZV/s0Njki37Pn/B4AGpRvYHIkQngntyUIrfVZ4Kxl+oZSaj9QBcj3CcL2UtaCBfwkObjJyMUjAemDECKn8qQPQikVAjQGNjtYNhwYDlCtWrW8CMc06WsryaWs7vVBxw/MDkEIr+b2BKGUKgrMB57XWl9Pv1xrPQWYAhAREaHdHY9ZtNZ8sHCn9bmUz3C/ZlWamR2CEF7NrQlCKRWAkRxmaa0XuPNYniouMZkHbMaOLlTAj59f6YK/n4wA527R54xih+EVw02ORAjv5M6rmBTwDbBfa/0/dx3Hk206+A/j5m6zPq9SughTn24nySGPPL/0eUD6IITIKXeeQdwDDAJ2K6XS6laP1VovzmIbnzF3/RGm/XEAgB5NqzGoXW1KFilkclT5y0ddPjI7BCG8mjuvYloH5Lufyimpmul/xvDjhiMAvDugBU2qlzU5qvxJmpaEyB25k9oFUlJT2Xn8MusOnGXJ9pOkaqOv/buR91GxVJDJ0eVfW09vBaSzWoickgSRA1prNsb8w56Tl5m/6ViG5Q+2COWxNjUpLjWVTPXyipcB6YMQIqckQdyBSzfiWbX7NMt2nOTU5VsABPj74e+nCClfjEHtatMopAwB/n4mRyoAPu32qdkhCOHVJEE44eL1eGb9dYilO4zmo9qVS/BM53o0CilLSPliZocnMiElNoTIHUkQWUhJ1SzefoJpqw6QmJxCl8ZV6dG0GjUqljA7NOGEDSc3AFK0T4ickgSRifjEZN6et52oIxcIDynDc90bUqV0EbPDEndg7KqxgPRBCJFTkiAcuBmfxH/nbGX/qSv8u1sDujephnHfn/AmX/X4yuwQhPBqkiDSuXorgbGztvD3hRuM7dOEe+tVMjskkUMy1KgQuSMJwsathCRembGZs1du8Wb/CBndzcutOb4GgHYh7UyORAjvJAnCxuTf93Di4k0mPNaMptXLmR2OyKVxq8cB0gchRE5JgrDYcOAcq/eeYUhkbUkOPmJa72lmhyCEV5MEgdEp/cmSPYSWL0a/1jXMDke4SPVS1c0OQQivJrf8AlNX7ufqrQRe6tWIAnIXtM9YeXQlK4+uNDsMIbxWvj+D2HHsIkt3nOThVtWpVUlugPMlE9ZOAKBD9Q4mRyKEd8rXCSIxOYWPf99NldJFGNSuttnhCBeb8eAMs0MQwqvl6wTx8+bjnL0Sy8SBLSgUIOND+5qqJaqaHYIQXi3fNrhfvhnP7HWHaFmrPI1DZUAfX7T08FKWHl5qdhhCeK18ewYx7Y8YkpJTGd6pntmhCDeZuG4iAF1qdjE5EiG8U75MEDFnrrJi5ykeblVdCvD5sDl955gdghBeLd8liFSt+WLpXkoXLcRj99YyOxzhRhWLVjQ7BCG8Wr7rg/hz92n2n77K4/fXIahQvsuP+cqimEUsillkdhhCeK189Q15Kz6Jb/44QJ3KJekQFmx2OMLNPtz4IQA96/Q0ORIhvFO+ShBfrzrA5RsJjOsXgZ+M7+Dz5vWbZ3YIQni1fJMgth+9yOLtJ3ioZSh1Kpc0OxyRB8oGyeXLQuRGvuiDiE1I5qPfdhFcughDImUQmfxiwf4FLNi/wOwwhPBa+eIMYurK/Zy/FseHQ1vJHdP5yOTNkwHoc3cfkyMRwjv5fILYevg8i7ef4OFW1alftbTZ4Yg8tPCRhWaHIIRX8+kEcT0ukf8t2sVd5YoyOFKK8eU3JQKlOq8QueHTfRBTVuznWmwiL/cOp2ABaVrKb+bumcvcPXPNDkMIr+WzZxA7jl1kxc5T9L+nhozzkE99se0LAPo36G9yJEJ4J59MEAlJxjgPlUsHMUDKaeRbiwcsNjsEIbyaTzYxzVp7iLNXYhnVraFctZSPBQUEERQQZHYYQngtn0sQB89c5aeNR+nUKJhwGechX5u5ayYzd800OwwhvJZPNTElJKXwwcKdlC5aiOEdZZyH/O7r7V8DMDBsoMmRCOGdfCZBaK3536JdnLh4k3cea06xwgFmhyRMtmLQCrNDEMKrubWJSSnVRSkVo5Q6rJQa467jJCSl8OGiXazee4Yn7q9D0xrl3HUo4UUC/AMI8JcfCkLklNvOIJRS/sBnQEfgFLBVKfWr1nrfnexHa02qNgb60VqTmqpJ0ZqUVM3ZK7Hs/vsyv0X9zdkrsQxsW4t+rWu44+UILzQ9ejoAQ8OHmhqHEN7KnU1MzYHDWuujAEqpOUBvINMEcfjcdXq9u8SaEFJTNdqJA9WtUpLnujWkSXXplBa3SYIQIneU1s58Bedgx0r1BbporZ+0PB8EtNBaj0y33nBguOVpA2CPWwJynbLARbODcILE6VoSp2tJnK5TR2tdzB07ducZhKMReTJkI631FGAKgFJqm9Y6wo0x5Zo3xAgSp6tJnK4lcbqOUmqbu/btzk7qU0BVm+fBwBk3Hk8IIYQLuTNBbAVqKaVClVIFgUeAX914PCGEEC7ktiYmrXWyUmoksAzwB6Zprfdms9kUd8XjQt4QI0icriZxupbE6Tpui9FtndRCCCG8m8/VYhJCCOEakiCEEEI45BEJIq9KcmRx/KpKqT+VUvuVUnuVUqMs899USp1WSkVbHt1stnnVEm+MUqpzXr0WpdRxpdRuSzzbLPNKK6VWKKUOWf4tZZmvlFKTLbHsUko1sdnPEMv6h5RSQ1wYXx2b9ytaKXVdKfW8J7yXSqlpSqnzSqk9NvNc9t4ppZpa/jaHLds6utQ7p3F+oJQ6YInlZ6VUScv8EKVUnM37+mV28WT2ml0Up8v+zsq4wGWzJc65yrjYxVVxzrWJ8bhSKtoy35T3U2X+HWTu51NbSliY9cDowD4CVAcKAjuBenkcQyWgiWW6GHAQqAe8CYx2sH49S5yFgFBL/P558VqA40DZdPPeB8ZYpscA71mmuwFLMO5JaQlstswvDRy1/FvKMl3KTX/bc8BdnvBeAm2BJsAed7x3wBaglWWbJUBXF8bZCShgmX7PJs4Q2/XS7cdhPJm9ZhfF6bK/M/Aj8Ihl+kvgaVfFmW75h8B/zXw/yfw7yNTPpyecQVhLcmitE4G0khx5Rmt9Vmu93TJ9A9gPVMlik97AHK11gtb6GHAY43WY9Vp6A99Zpr8DHrCZ/702bAJKKqUqAZ2BFVrry1rrK8AKoIsb4moPHNFa/51N7HnyXmqt1wKXHRw/1++dZVlxrfVGbfxv/N5mX7mOU2u9XGudbHm6CeO+okxlE09mrznXcWbhjv7Oll+39wPz3Bmn5Tj9gNlZ7cPd72cW30Gmfj49IUFUAU7aPD9F1l/ObqWUCgEaA5sts0ZaTuGm2Zw6ZhZzXrwWDSxXSkUpo0wJQAWt9VkwPmhAeQ+IE4x7X2z/43naewmue++qWKbdHS/AExi/ANOEKqV2KKXWKKXutczLKp7MXrOruOLvXAa4apMU3fV+3gv8o7U+ZDPP1Pcz3XeQqZ9PT0gQTpXkyAtKqaLAfOB5rfV14AugBhAOnMU4FYXMY86L13KP1roJ0BV4VinVNot1TYvT0l7cC/jJMssT38us3GlceRKvUuo1IBmYZZl1FqimtW4MvAj8oJQqnlfxOOCqv3Nexf8o9j9iTH0/HXwHZbpqJvG49P30hAThESU5lFIBGH+YWVrrBQBa63+01ila61RgKsbpMGQes9tfi9b6jOXf88DPlpj+sZxCpp0Knzc7TowEtl1r/Y8lXo97Ly1c9d6dwr7Zx+XxWjocewADLM0EWJpsLlmmozDa82tnE09mrznXXPh3vojRbFIg3XyXsey7DzDXJn7T3k9H30FZ7DtvPp932pni6gfG3dxHMTqu0jqp6udxDAqjTe6jdPMr2Uy/gNGGClAf+w63oxidbW59LUARoJjN9AaMvoMPsO/Iet8y3R37jqwt+nZH1jGMTqxSlunSLn5P5wCPe9p7SbpOSFe+dxjlZVpyuxOwmwvj7IJRKr9cuvXKAf6W6erA6eziyew1uyhOl/2dMc4+bTupn3FVnDbv6RpPeD/J/DvI1M+ny74QcvPA6JE/iJGtXzPh+G0wTrd2AdGWRzdgBrDbMv/XdB/+1yzxxmBzNYA7X4vlA7vT8tibtn+M9tpVwCHLv2kfCIUxaNMRy+uIsNnXExgdhYex+SJ3UZxBwCWghM08099LjKaEs0ASxi+qf7nyvQMiMMrVHwE+xVKpwEVxHsZoW077fH5pWfchy2dhJ7Ad6JldPJm9ZhfF6bK/s+XzvsXy2n8CCrkqTsv86cBT6dY15f0k8+8gUz+fUmpDCCGEQ57QByGEEMIDSYIQQgjhkCQIIYQQDkmCEEII4ZAkCCGEEA5JghBeSSm1Winl9sHklVLPWSpszko3P0IpNdkyHamUau3CY4YopR5zdCwh8pLbhhwVwlMppQro2zV+svMMxjX7x2xnaq23AdssTyOBmxg3LroihhDgMeAHB8cSIs/IGYRwG8sv4f1KqamWGvfLlVKFLcusZwBKqbJKqeOW6aFKqV+UUouUUseUUiOVUi9aiqdtUkqVtjnEQKXUBqXUHqVUc8v2RSxF4rZatults9+flFKLgOUOYn3Rsp89SqnnLfO+xLhZ61el1Avp1o9USv1mKaz2FPCCMsYPuFcpVU4pNd8Sw1al1D2Wbd5USk1RSi0Hvre8P38ppbZbHmlnIROBey37eyHtWJZ9lLa8P7ss70eYzb6nWd7Xo0qp52zej9+VUjstr61/7v6qIl/J7d2p8pBHZg+MX8LJQLjl+Y/AQMv0aix3fwJlgeOW6aEYd4AWwyh7cA3L3a7A/8MoYpa2/VTLdFssZRSAd2yOURLjDt0ilv2ewsFdrkBTjLtRiwBFMe6kbWxZdpx0429Y5kcCv1mm38RmDASMX/5tLNPVgP0260UBhS3Pg4BAy3QtYFv6fTs41ifAOMv0/UC0zb43YJSyKItxJ3sAxp3BU232VSL9a5GHPDJ7SBOTcLdjWutoy3QURtLIzp/aqIl/Qyl1DVhkmb8bCLNZbzYY9f6VUsWVMcpaJ6CXUmq0ZZ1AjC9psNTJd3C8NsDPWutbAEqpBRhloHc48wId6ADUU7cH7CqulCpmmf5Vax1nmQ4APlVKhQMpGEXhstMG40sfrfUfSqkySqkSlmW/a60TgASl1HmgAsZ7Nkkp9R5Gkvkrh69J5EOSIIS7JdhMpwCFLdPJ3G7iDMxim1Sb56nYf2bT14lJK2v8kNY6xnaBUqoFcCuTGHM0NGgW/IBWNokgLQbSxfAC8A/QyLJNvBP7zqpsc/r3uoDW+qBSqilGXZ93lVLLtdZvO/UqRL4nfRDCLMcxmnYA+uZwH/0BlFJtgGta62vAMuDfSlnHC27sxH7WAg8opYKUUkWAB4E7+aV9A6NJLM1yYGTaE8sZgiMlgLPaKI09CKO6qaP9pY91gGW/kcBFncW4AUqpykCs1nomMAlj6E0hnCIJQphlEvC0UmoDRpt5TlyxbP8lRiVRgPEYTTe7lDFI/fjsdqKNoR6nY1QO3Qx8rbW+k+alRcCDaZ3UwHNAhKUjeR9GJ7YjnwNDlFKbMJqX0s4udgHJlo7lF9Jt82bavjE6s4dkE1tDYItSKhqjmuqEO3hdIp+Taq5CCCEckjMIIYQQDkmCEEII4ZAkCCGEEA5JghBCCOGQJAghhBAOSYIQQgjhkCQIIYQQDv1/ZNiXijrZtRcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.axhline(mi,label='ground truth',linestyle='--',color='red')\n",
    "for i in range(rep):\n",
    "    plt.plot(mi_list[i,:],color='steelblue')\n",
    "    for t in range(mi_list[i].shape[0]):\n",
    "        if (mi_list[0,t]>1*mi):\n",
    "            plt.axvline(t,label='100% exceeded',linestyle=':',color='green')\n",
    "            break\n",
    "plt.xlim((0,mi_list[0].shape[0]))\n",
    "plt.ylim((0,mi*2))\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"MI estimate\")\n",
    "plt.legend()\n",
    "plt.savefig(fig_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
